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SUMMARY 


This  report  presents  an  analysis  of  a DC  power  supply  consisting 
of  a superconducting  alternator,  a rectifier  bridge,  and  an  LC  output 
filter.  The  main  purpose  of  this  researcli  was  to  determine  if  changes 
in  the  size  of  the  alternator  inductances  would  allow  the  use  of  a smaller 
filter.  To  perform  this  study  it  was  necessary  to  examine  the  be- 
havior of  the  filter  and  to  determine  how  its  operation  was  affected  by 
the  alternator  parameters. 

Basically,  the  filter  performs  two  functions: 

1.  It  attenuates  the  output  ripple  voltage. 

2.  It  limits  the  initial  fault  current  when  a short  circuit  occurs 
at  the  load. 

Both  of  these  functions  also  depend  upon  the  values  of  the  alternator 
inductances . 

Since  the  first  function  refers  to  the  stead'  state  behavior,  it 
was  necessary  to  develop  a model  for  this  operating  mode.  This  was  done 
first  for  a system  with  an  uncontrolled  rectifier  bridge  and  then  these 
results  were  extended  to  a controlled  rectifier  bridge  system.  The 
second  function  is  a transient  phenomenon,  so  it  was  also  necessary  to 
develop  a second  model  to  describe  the  transient  behavior. 

Once  the  system  models  were  complete,  a study  was  performed  where 
the  unfiltered  ripple  voltage  w.is  calculated  for  various  values  of 
the  alternator  inductances.  It  was  found  that  under  certain  conditions 
the  ripple  voltage  can  be  decreased  by  increasing  the  armature  self 


iv 


inductance  (L  ).  A program  was  then  written  which  calculated  the  weight 
of  the  LC  filter  that  was  require-d  for  a given  set  of  specifications  and 
alternator  parameters.  This  program  indicated  that  an  increase  in  L 

3 

could  decrease  the  required  filter  weight  by  as  much  as  22%. 

Other  investigations  included  a sensitivity  analysis  of  the  alternator 
inductances  and  the  design  and  testing  of  a phase  controlled  voltage 
regulator  with  current  overload  protection. 
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1.  ALTERNATOR  WITH  UNCONTROLLED  RECTIFIER  BRIDGE 


1.1  Introduction 

Recent  advancements  in  superconducting  alternators  have  created  a 
strong  interest  in  using  these  machines  for  airborne  electric  power  sup- 
plies. The  predominant  advantage  of  this  power  source  is  its  relatively 
low  weight  for  applications  requiring  multi-megawatt  outputs  at  several 
kV.  This  low  weight  characteristic  occurs  because  of  two  factors: 

1.  Even  with  the  required  cryogenic  equipment,  the  super- 
conducting alternator  system  weighs  much  less  than 

a conventional  alternator. 

2.  The  higher  armature  voltages  of  the  superconducting 
machine  may  eliminate  the  need  for  heavy  output 
inverters  and  transformers. 

These  attributes  are  discussed  in  further  detail  in  such  references  as 
[1]  - [12],  and  a very  recent  example  of  such  a machine  is  described  by 
McCabria,  et  al.  in  [13,14].  This  particular  machine  develops  lOMVA  at 
5 kV  and  weighs  approximately  1,000  pounds  (alternator  weight  only). 

This  same  reference  also  includes  projected  estimates  for  a 25MVA  machine 
weighing  between  1882  and  2160  pounds,  depending  on  rated  output  voltage 
(again,  these  figures  only  include  the  weight  of  the  alternator). 

The  potential  advantages  of  superconducting  alternators  have  prompted 
extensive  research  in  this  area,  mo: ; of  which  has  concentrated  on  ac 
loads  (again  see  [1]  - [14]).  Applications  for  these  machines  also  exist 
in  high  power  dc  systems  however,  where  the  alternator  is  connected  to  a 
rectifier  bridge  followed  by  a large  filter  choke.  This  mode  of  operation 
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has  been  studied  in  detail  for  convent io''al  alternators,  (see  [15]  through 
[21]),  but  until  now  no  such  analysis  has  been  presented  for  the  super- 
conducting machine. 

One  of  the  more  rigorous  analyses  of  conventional  rectified  alter- 
nators is  that  presented  by  Franklin  [17,18]  for  salient  pole  machines. 

By  assuming  constant  flux  linkages  for  the  rotor  windings,  this  study 
derives  a set  of  nonlinear  equations  in  terms  of  the  electrical  variables 
of  interest.  Certain  approximations  then  lead  to  a linearization 
involving  a constant  K factor,  and  an  explicit  solution  is  obtained.  The 
advantages  of  this  approach  are  readily  apparent  since  it  provides  a 
closed  form  expression  for  each  of  the  variables,  once  the  proper  K factor 
has  been  found.  The  determination  of  K is  somewhat  distracting  however, 
since  it  is  load  dependent  and  requires  the  use  of  numerical  methods.  In 
the  following  section  it  will  be  shown  that  this  K factor  can  actually  be 
eliminated  from  the  final  solution  if  a Newton-Raphson  algorithm  is  used. 
This  new  approach  appears  to  have  certain  advantages  since  it  is  somewhat 
less  complicated  and  does  not  depend  on  any  linearization  factors. 

The  essence  of  the  work  presented  here  is : 

1.  Franklin's  basic  analysis  methods  are  extended  to  the 
superconducting  machine. 

2.  The  dependence  on  the  previously  mentioned  K factor  is 
eliminated.  As  stated  above,  this  is  accomplished 

by  using  a Newton-Raphson  algorithm  where  a K=1 
is  used  only  to  find  a starting  point. 

3.  A numerical  example  predicting  the  rectified  charac- 
teristics of  the  machine  described  by  McCabria,  et  al. 
in  [13,14]  is  included. 
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The  overall  intent  is  to  provide  an  analytical  model  of  the  steady  state 
behavior  of  the  superconducting  alternator  with  a rectified  output. 

This  analysis  is  regarded  as  a preliminary  step  to  the  eventual  design 
and  testing  of  these  machines  for  D.C.  loads. 

1.2  Steady  State  Alternator-Rectifier  Model 

The  armature  of  the  supei'conducting  machine  is  assumed  to  be  Y 
connected  as  indicated  for  the  basic  two  pole  machine  in  Figure  1.  The 
d and  q windings  shown  in  this  figure  are  equivalent  windings  that 
account  for  the  effect  of  the  cylindrical  damper  shield  located  between 
the  rotor  and  the  stator  (see  [5],  [13]  or  [14]  for  example).  Output 
voltage  and  current  waveforms  are  shown  in  Figure  2,  where  the  indicated 
6 corresponds  to  Figure  1.  Formulation  of  this  problem  proceeds  in  much 
the  same  manner  as  in  [17,18],  but  there  are  some  important  differences 
in  the  machine  parameters.  It  also  should  be  noted  that  the  method  of 
solution  is  quite  different  from  these  earlier  refer<2nces,  and  certain 
equations  are  employed  in  a different  manner. 

The  following  approximations  are  utilized: 

1.  All  winding  and  diode  resistances  are  quite  small  and  can  be 
ignored. 

2.  All  diode  voltage  drops  are  negligible. 

3.  The  load  inductance,  L^,  is  sufficiently  large  to  maintain  a 

constant  Ij^,  i.c.,  the  effect  of  load  current  variations  is 
ignored. 

4.  Each  armature  winding  is  assumed  to  have  a perfect  sinu- 
soidal distribution  about  the  stator. 


3. 


Figure  1.  Equivalent  Circuit  for  the  Superconducting 

Alternator  with  Uncontrolled  Rectifier  Bridge. 


Figure  2.  Output  Voltage  and  Armature  Currents 
with  Uncontrolled  Rectifier  Bridge. 
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5. 


The  effect  of  the  damper  shield  is  modelled  by  equivalent 
direct  axis  and  quadrature  axis  windings  on  the  rotor 
(d  and  q). 

6.  The  rotor  speed  is  assumed  constant. 

Since  the  superconducting  alternator  is  an  air  core  machine  there  are  no 
saturation  or  saliency  effects. 

Although  the  line  to  line  voltages  in  Figure  2 are  shown  as  per- 
fect sinusoids,  it  should  be  noted  that  they  are  actually  distorted 
somewhat.  Thus  commutation  actually  starts  at  some  B > 90°  and  not  at 
9=6=  90°  as  indicated  in  the  figure.  Another  interesting  characteristic 
is  the  fact  that  the  stator  MMF  is  constant  in  magnitude  and  direction 
during  the  conduction  interval  and  abruptly  shifts  to  a new  direction 
during  the  commutation  interval.  This  phenomenon  is  termed  "MMF  jump" 
as  is  described  further  in  such  references  as  Stepina  [16]  and  Franklin 
[17]. 

The  waveforms  shown  in  Figure  2 indicate  that  p < tt/3.  Franklin 
points  out  that  it  is  also  possible  to  reach  a mode  where  = "J"  Although 
a large  number  of  simulations  were  conducted  in  this  present  study,  the 
p = 7t/3  mode  was  never  reached  for  the  machine  used  in  the  numerical  exam- 
ple. The  conclusion  drawn  was  that  this  appeared  to  be  an  unlikely  operat- 
ing mode  for  this  application,  so  it  was  not  included  in  the  analysis. 

1.3  Steady  State  Equations 

The  primary  goal  of  this  section  is  to  derive  five  equations  that 
are  expressed  in  terms  of  the  following  variables; 
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P = angle  at  which  commutation  starts 

p = commutation  angle 

Ij-  = average  field  current 

W = variable  defined  by  equation  (20.) 

V = variable  defined  by  equation  (21.) 


These  equations  turn  out  to  be  nonlinear  with  respect  to  3 and  w, 
but  they  can  be  solved  by  some  numerical  method  such  as  the  Newton- 
Raphson  algorithm.  Once  these  variables  have  been  found  it  is  possible 
to  determine  the  time  dependent  expressions  for  the  output  voltage  and 
the  current  in  each  winding. 

It  is  assumed  that  the  field  current  consists  of  the  constant  com- 
ponent, I^,  and  a time  varying  component,  i^. 


^f(tot)  ' ^f  ^f 


(1.) 


The  winding  currents  during  conduction  and  commutation  are  indi- 


cated as  follows. 

Conduction  (Interval  1-2  in  Figure  2),  6 + u - tt/3  ^ 9 < 3 


— 
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Commutation  (Interval  2-3  in  Figure  2),6<0<3+p 


'-'l  * ‘k> 


(If  + if) 


The  flux  linkages  are  given  by  the  following  expression, 

r n 


= [L]  i 
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Figure  2 indicates  that. 


V = V = -0)  — [X  - X ],  6+W  - ir/3  < w t < B+u 
o ah  d0  a b 


The  average  voltage,  V , at  the  output  of  the  rectifier  bridge  is, 


= I ^^023 

B+u-it/3 


where  ^nd  represent  the  output  voltage  functions  over  1- 


2-3  respectively. 


V = - ^ 
L It 


; 


B 


B+m-h/3 


dX  dX^  „ dX  dX, 

[jT  - Car  - 


Figure  2 indicates 


V = V 0 0 = 6 

012  023 


V = - ^ 

L It 


^b^23 


- (X 


(6+p) 


^^12 


(B+Vi-tt/S) 
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(6.) 

(7.) 

(8.) 

2 and 

d0  (9. ) 

(10.  ) 
(11.) 


(2.)  and  (3.)  become  the  same  when  ij^  = 0,  therefore  using  (3.)  and  (5.) 


X = ' I M sin  (0+ti/6)  + i M cos0  + i L. 
q L d k d q d 


where  i,  = 0 over  the  conduction  interval, 
k 


X is  assumed  to  be  constant  over  1-3.  At  6 = 6+u-n/3  we  have  i = i , 
q q qo’ 


X (6+vi-Tt/3)  = - /3  I M sin  (B+y-ir/6)  + i L, 
q L d qo  d 

i = - /s  I,K  [sin(  B+u-ir/6)-sin(  0+tt/6)  3-/3  i,  K cos0+i 
q L q k q qo 

where  K = M,/L, 

q d d 

Using  a similar  procedure,  expressions  for  X^  and  X^  may  be  determined. 

The  two  simultaneous  equations  for  X^  and  X^  may  then  be  solved  for  i^ 

and  i , , 
d 

i,  = I,  K , [cos  (B+u-w/6)  - cos  (0+71/6)3  - /s  i,  K , sin0  + i , 
d L d k d do 

i_  = /T  [cos  (B+u-tt/6)  - cos  (0+ir/6)3  - i,  K sinO  + i^ 
t L f k t to 


where 


"/■d  - "d"fd 
‘■/'d  - 


Vf  - "f»fd 
•■f^d  - <"fd>' 


As  pointed  out  by  Shilling  [153,  the  rotor  currents  are  periodic  with 
respect  to  the  6th  harmonic;  therefore. 


5 i^dB  = / i,d6  = / 1 do  = 0 

f d q 


B+u-ii/3 


B+u-tt/3 


B+w-tr/3 


Since  ij^  = 0 for  6+m  - y < 0 < B,  we  may  define  the  following  constants 


W i.  i.  sin0  d6 

B 


11. 


(12.  ) 

(13.  ) 
(14.  ) 
(15.) 

(16. ) 
(17.) 

(18.  ) 

(19.) 

(20.) 


V 


i,  cos6d6 
k 


(21.  ) 


- 

6 K 

Integrating  (14.),  (16.)  and  (17.)  over  6+p-  j - 9 - B+u  and  solving 

for  i , i,  and  i,  produces, 
qo’  do  fo  ^ 

i = /b  I K [sin  (B+u-ti/6)  - — sin  (B+w)]  + K V 
qo  L q TT  ^ <1 

i,  = -/b  I,K  [cos  (B+y-iT/6)  - — cos  (B+u)]  + K W 
do  L d ■rr  it  d 

*^f 

^fo  ■ ^do 

d 

Therefore,  substituting  into  (14.),  (16.)  and  (17.), 


i = /B  K {—  [V  - I, sin  (B+u)]  + [I  sin  (0+tt/6)  - i cos  (0)]} 
q q n L L k 

i . = /b  K {—  [W  + I cos  (B+w)]  - [I, cos  (0+11/6)  + i,  sin  (0)]} 

cl  Q TT  Ij  L>  K 

^f  ^d 

d 

Drom  ( 4 . ) , we  have , 

X = I,  (L  + M ) + [(I^  + i^)  + i.M,]  cos0  - i M sin0 

aLaa  fffdd  qd 

X,  = - I-  (L  + M ) + i,  (L  + M ) + [(I.  + i.)  + i,M,] 

D Laa  kaa  fffdd 

cos  (0  - 1^)  - i M sin  (0  - |^) 

3 ^ d 3 

Utt 

X = - i,  (L  + M ) + [(I^  + i ) M_  + i.M.]  cos  (0  - — ) 
c Kaa  fffdd  B 

- i M sin  (0  - ^) 

q d 3 

Using  the  results  of  (25.)  - (27.), 

(I^  + i^)  M-  + i.M,  = I.M.  + — M W - /b  I,  M [cos  (0+n/6) 
r rr  dd  rr  tto  Lo 

- — cos  (B+u)]  - i/T  i,  M sin0 

TT  K O 


(22.  ) 
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(29.  ) 
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i M 
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^ M V + /3  I M [sin(0+Ti/6)  - - sin  ( 6+U ) ] 
™ oo  L cx) 

- i,  M cos  6 

k oo 


where,  M = (K^M^  + K.M,) 
o f f d d 


‘■d'-f  - <"fd'' 


M = K 
oo  q d 


(11.)  can  be  used  to  find  V^,  by  noting  that 


i^=I^@  0=6+p 


i , = 0 @ e = B +U-1T/3, 
k 


Vl  = ^ II  *o  + |-  cos  (20+2U+1T/3)  + I^M^in(6+u) 


It  a.  sin  (2B+2w)  + — [M  Wsin(6+w)  + M Vcos  (B+u)]) 
2tt  L d It  o oo 

where,  A^  = (M  + M ),  A,  = (M  - M ) 
f o oo  d o oo 

A = 1 (L  + M ) - A- 

o 3 a a f 


The  current  ij^  exists  only  during  the  commutation  period  where  the 
"b"  and  "c"  phases  are  shorted  together,  i.e. , 


= 0,  e - 0 - 8+p 


(X,  - A ) = constant,  6 - 0 - 6+u 

D c 

For  0 = 3,  i,  = 0,  therefore  setting  (A,  - A ).  = (A  - ^ )q  one  ob- 

K D C O D C p 

tains , 


. _ 1 

^k  Ta  + A ,cos( 20)] 
o u 


2MfI  g 

(sinB-sin0)  + — M W(sin6-sin0) 


(32.  ) 
(33.  ) 

(34.) 


(35.) 

(36.) 


13. 


[sin(26+lJ)  - sin  (6+m+0)] 


F,  31  A 

+ — M V (cos6-cos9)  + 

It  OO 


[sinp  + sin  (0-0-p)]  - 


IjA,  [cos  (26-71/3) 
L d 


- cos 


(20-71/3)] 


(37.  ) 


Again  utilizing. 


be 


0,  6-0-6+y 


we  have  for  i,  = 0 @ 0 
k 


- ''c^a 


= e, 

- ‘c>6  ' “ 


(38.  ) 


which  leads  to, 

1 3 

-A  sin  (26  - tt/3)  = — — cos  6 + (M  W cos  6 - M Vsin6' 

d /Jh  ^ '’h  ° °° 


+ [A^cos  (26+y)  + A^  cos  u]  (39.) 

Utilizing  ij^  = @ 0 = 6+u  in  (37.)  leads  to, 

I 

+ 2A,  cos  (26+lJ)  cos  (m+ti/3)  = - — M [cos(6+  ^)  sin  (■^)] 

+ — ^ tM  W(sin6-sin( 6+P  ) ) + M V (cos6-  cos(6+u))] 

17 1 o oo 

- f [2A^  cos(26+  -|^)sin  (j)  + A^sin(u)]  (40.) 

One  could  substitute  (37.)  into  (20.)  and  (21.)  and  integrate  to  find 
two  more  equations,  which  along  with  (35.),  (39.)  and  (40.)  would 
yield  five  nonlinear  equations  for  the  five  unknowns,  I^,  6,  u,  V and  W. 
This  process  is  simplified  considerably  by  use  of  the  following  approxi- 
mation , 


A,  = 0,  (i.e.jM  »M  ) 

d o oo 


(41.  ) 


14. 


Equations  (35.),  (37.),  (39.)  and  (40.)  indicate  that  always  appears 
in  conjunction  with  or  A^.  Therefore,  (41.)  is  acceptable  if  A^  is 
small  in  comparison  to  A^  and  A^. 

The  superconducting  alternator  considered  in  this  study  (the  same 
machine  described  by  McCabria,  et  al. , in  [13,14]  has  the  following 
parameters^ : 


^f 

= 1.2 

H. 

= 7.9  X 10  ^H. 

"d 

= 8.2  X 

10"®H. 

= 1.9  X lO” 

ru 

L 

= 3.0  X 

10"‘*H. 

= 3.8  X lO”^ 

a 

d 

M = 1.5  X 10"‘^H. 

a 

A_.  = 0.01  X lO"*^  H.,  A^  = 3,5  x lO"**  H.  , A = 2.5  x 10“*^  H.  (43.) 

d to 

Therefore  the  approximation  given  by  (41.)  appears  to  be  acceptable, 

at  least  for  this  particular  example. 

Using  (41.),  equations  (20.),  (21.),  (35.),  (39.)  and  (40.)  reduce 


to, 

0 = A W + A (cos  (6+p)  - cos  6)  + ^ ( 2u-sin( 26+2p )+sin26 ) 
o 4 

+ j (sin^(6+y)-sin^e)  (44.) 

0 = A V + A(sinB-sin( B+y ) ) + ^ (sin^( 6+y  )-sin^B ) 
o 2 

+ ^ (2y+sin(26+2y)  - sin(2B))  (45.) 

0 = -V,  + — [-  I,A  + /b  I_M.sin  (B+y) 

L Tf  ‘-4  L o f f 

9M 

+ — ^ (W  sin  (B+y)  i V cos  (y+y))]  (46.) 

7T 


These  parameters  were  supplied  by  H,  Southall  of  the  U.S.  Air  Force 
Aero  Propulsion  Laboratory. 
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0 


cos  6 + 


(47.) 


3M 


/a  I, 


°°  (W  cos  B-VsinB)+— M cos  4 


IT  oo 


4 I 


0 = -A 


° ^ I 


^ 6 M 6 M 

f /•  o 4 . . , U ^ oo  . OO 

— M..  cos  (B  + tt)  sin  (— ) - sin  u + — = — . 

r 2 2 T(  TTl 


[W  (sin  B - sin  (6+u))  + V (cos  B - cos  (B+u))] 


(48.  ) 


where 


2 6 M 

A = — sin  6 + — (W  sin  6 t V cos  B - I,  sin  u) 

^ f f 1T  L 


(49.  ) 


B = — + - M (W  + I,  cos  (B  + u)) 

^ f f TJ  oo  L 


(50.  ) 


6 M 


C = 


oo 


(V  - I,  sin  (B  + u) ) 
TT  L 


(51.  ) 


(44.)  - (48.)  provide  five  nonlinear  equations  which  are  functions  of 
the  variables  B,  4,  V and  W.  Actually,  these  equations  are  linear 
with  respect  to  1^;,  V and  W so  it  would  be  possible  to  eliminate  these 
variables  and  have  a set  of  two  equations  which  are  functions  of  B and 
4.  The  equations  involved  in  this  reduction  are  quite  cumbersome  how- 
ever, so  one  might  as  well  work  directly  with  (44.)  - (48.). 


1.4  Solutii'n  for  I^,  B,  y , V and  W 

Rewriting  the  variables  and  equations  in  matrix  form, 

6 

4 

X ^ If 

V 


(52.  ) 


W 


f (x)  = R.H.S.  of 


(44.  ) 


(46.) 
(47.  ) 


(48.  ) 


3f  ( x) 


Jacobian  matrix  = F (x)  = 


(52.)  - (54.)  can  be  used  to  form  the  standard  Newton-Raphson  equation, 

F (x  ) (x  - x ) = f (x)  - f (x  ) ([ 

— o — — O — _ — — O 

where  ^ is  some  initial  starting  point  which  must  be  reasonably  close 
to  the  desired  x.  In  this  case  f (x)  = 0,  so  we  have. 


F (x  ) (x  - X ) = - f (x  ) 

It  remains  to  find  a satisfactory  value  for  This  is  accomplished 

by  making  use  of  the  linearization  suggested  by  Fr.nklin^, 


L = (e  - 6)  — 

k p 


Substituting  (57.)  into  (20.)  and  (21.)  gives  the  result. 


W = -Ij^  cos  ( 8+p ) + 


2 I^  sin  (w/2)  cos  (6+p/2) 


Actually,  Franklin  uses  a "K  factor"  as  mentioned  in  the  Introduction, 

II 

where  0,5  -K  - 0.9.  This  produces  the  approximation,  ij^  ==  K(6-6)  — . 

Since  (57.)  is  only  used  to  find  a starting  point  for  the  Newton-Raphson 
algorithm,  K is  not  critical,  and  K = 1.0  is  used. 


17. 


V 


I sin  ( 6+U  ) 
L» 


(59.  ) 


2 I sin  (Vi/2)  cos  (P+U/2) 

Ij 


After  substituting  (57.)  - (59.)  into  (35.),  (39.)  and  (40.)  and  per- 
forming some  rather  laborious  calculations,  one  obtains. 


-1 

U = cos 


4ttV,  - 9ul  A 
L L o 

4tiV,  + 9wl  A 
L L o 


e 


tan 


-1 


2 

ttA  12  M (sin  ^)  (1  - cos  u) 

o oo 2 

12  M sin  (u/2)  cos  (^)  (1  - cos  u) 
oo  2 


(60.  ) 


(61.  ) 


I 


f 


1 

/3  sin  (B+p) 


3 

4 


18  I,  M 
L oo 

np 


(sin  -) 


(62.  ) 


(58.)  - (62.)  provide  a value  of  which  produces  convergence  within 
three  or  four  iterations,  depending  on  the  convergence  tolerance. 


1.5  Numerical  Results 

This  example  is  based  on  the  same  4 pole,  400  Hz,  10  MVA/51<V  super- 
conducting alternator  described  by  McCabria,  et  al.  in  [13,14].  When 
operating  into  a bridge  rectifier  at  full  load  with  a large  filter 
inductance,  this  system  will  have  the  output  values  (see  [24]), 

V = 6760  V.  dc 

^ (6- 

I,  = 1420  A.  dc 

Li 

The  following  results  assume  that  the  system  is  operating  with  a closed 
loop  controller,  i.e.,  1^.  is  varied  to  maintain  a constant  V^. 
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The  inductance  parameters  for  this  machine  are  indicated  in  (42.). 

All  data  is  presented  in  terms  of  the  actual  magnitudes,  but  a per  unit 
system  would  serve  just  as  well. 

If  vs. 

Figure  3.  shows  vs.  where  I^  is  varied  to  maintain  a con- 

stant V^.  As  seen  from  the  curve,  the  variation  in  is  approximately 
linear  up  to  200%  of  the  full  load  value  of  Ij^. 

B and  p vs. 

Figures  4.  and  5.  indicate  the  variation  in  6 and  p respectively 
with  respect  to  1^^.  Figure  5.  indicates  that  p remains  less  than  60° 
as  mentioned  earlier. 

Harmonic  Content  of  v vs.  I.  ; 

o L 

One  of  the  more  important  problems  in  this  type  of  power  system  is 

the  weight  of  the  output  filter  (L^  and  as  show  in  Figure  1.).  In 

order  to  minimize  the  combined  weight  of  these  components  it  is  necessary 

to  know  the  harmonic  content  of  v under  all  load  conditions.  The 

o 

output  voltage  is  obtained  from  (7.)  by  substitution: 

Conduction  period  (ij^  = 0),B  + p-tt/3<0^B, 

v-,_  = u)  (/s  1^  M + — M W)  sin  (0  + ti/6)  t ~ M V cos  (0+ii/6) 

012  f f TI  OO  IT  oo 

9w  I M 

+ ^ sin  (6+it/6-B-p)  (65.) 

IT 

Commutation  period  (ij^^O),  B<0<  B+M, 
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0.00  S.OO  10.00  15.00  20.00  25.00 

(>10'  J 

II  (amps) 


Figure  3.  versus  with  Uncontrolled  Rectifier 
Bridge . 
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o 


0.00 
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10.00  15.00  20.00  25.00 

(>10^  ) 

II  (amps) 


} 


Figure  4. 


6 versus  I with  Uncontrolled 
Kectifi"-'  Bridge. 
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.00  1.00  2.00 


Figure  5,  p versus  I with  Uncontrolled 

L 

Rectifier  Bridge. 
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r 


f 


,,  3M  2M^I-  6M  W 

^023  = ^^012  ^ a;  ^ ^ -^>  ° 


6M  V 61  M 

OO  -a  . Loo  /a  o 
Sin6  + COS  (0-B-u)] 


(66.  ) 


The  Fourier  series  for  is  given  by. 


= V + a cos(no)t)  + I b sin(ncjt)  n = 6,12,18,, 

o b - n r-  n 

n=6  n=6 


(67.) 


where  V = — v dO  + — V-„_  d 6 

L Ti  012  71  023 

6+ u- 71/3  8 


(68.  ) 


6 6 h 

a = — / V cos(n8)  d0  + — / v cos(ri0)  d0 

" " B+M-1T/3  " 6 


(69.  ) 


6 8 B+  p 

b = - / V sin(n0)  d0  + - / v.„_  sin(n0)  d0 

n 77  012  77  023 

B+u-77/3  6 


(70.) 


The  peak  value  of  the  nth  harmonic  is  given  by, 

. 2 2/^^^ 

c = (a  + b ) 
n n n 


(71.  ) 


Although  tedious,  the  evaluation  of  these  coefficients  is  straightforward 
and  will  not  be  included  here. 

The  peak  value  of  the  first  three  ac  harmonics  of  v vs.  I,  are 

o L 

indicated  in  Figure  6.  As  would  be  expected,  the  ac  content  is  dominated 
by  the  6th  harmonic. 


Variation  in  i,,  i,,  and  i : 
f d q 

As  pointed  out  in  various  references  ([8]  through  fll],  for  example), 
one  of  the  more  critical  problems  in  a superconducting  machine  is  heating 
due  to  induced  currents  in  the  field  winding  since  this  may  cause  the 
field  to  depart  from  the  superconducting  mode.  The  thermal  analysis 
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09 / 00  O 


Figure  6.  First  Three  Harmonics  of  v versus 

o 

I with  Uncontrolled  Rectifier  Bridge 


of  this  winding  is  bevond  the  scot>e  of  this  report,  but  it  is  certainly 


related  to  the  variations  in  i,.  The  instantaneous  value  of  i,  is 

f t 

qiven  bv  (27.)  and  plotted  as  a function  of  6 in  Fiqure  7.  The  rms 
value  of  ij  was  obtained  bv  numerical  integration  and  is  shown  as  a 
function  of  I in  Figure  8. 

Since  ladDoratorv  data  on  am  actual  rectified  superconducting  alter- 
nator was  not  available  at  the  time  of  this  report,  it  is  difficult  to 
say  how  these  calculated  variables  will  eventually  compare  with  experi- 
mental results.  However,  an  analytical  review  of  the  i^  calculations 
does  indicate  that  the  ij's  shown  in  Figures  7 and  8 may  be  higher  than 
the  actual  values.  It  is  believed  that  the  reason  for  this  is  that  the 
inductance  matrix  in  (5.)  does  not  fully  account  for  the  high  frequency 
attenuation  provided  by  the  damper  shield.  This  effect  is  still  under 
investigation. 
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Rectifier  Bridge.  (See  text  for  discussion 
of  i_  calculations.) 
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Figure  8.  RMS  value  of  i versus  I with  uncontrolled 

r L 

Rectifier  Bridge.  (See  text  for  discussion  of 
i^  calculations.) 
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Armature  Current  vs  9: 


The  variation  in  i^  during  conduction  and  commutation  is  shown  in 
Figure  9.  Note  that  the  commutation  value  is  simply  - i^^  as  given 
by  (37.). 


I 


28. 


I n 

7.00 


I I I 1 

8.50  10.00  11. SO  13.00  1 4 . SC 

(>10‘  ) 

0 (degrees) 


Figure  9. 
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2.  ALTERNATOR  WITH  CONTROLLED  RECTIFIER  BRIDGE 


2.1  Introduction 

It  is  probable  that  most  rectified  alternator  applications  will 
require  some  type  of  closed  loop  voltage  regulation.  A good  indication 
of  this  is  provided  by  McCabria,  et.  al.  [13],  which  describes  a 3 
phase  10  MVA/5kV  machine  that  has  an  open  loop  voltage  regulation  of 
26.5%.  The  previous  analysis  is  applicable  for  the  most  obvious  means 
of  control,  which  is  to  vary  the  current  in  the  superconducting  field 
winding.  This  section  considers  another  possibility,  which  is  to 
control  the  turn-on  angle  of  the  output  rectifier  bridge  ( in  this  case 
composed  of  thyristors). 

Although  phase-control  is  widely  used  when  an  A.C.  bus  is  the 
source,  this  method  is  usually  not  employed  with  a dedicated  rectified 
alternator.  In  the  case  of  a conventional  alternator  it  is  much  better 
to  regulate  the  voltage  by  means  of  field  control  since  this  technique 
is  relatively  simple  and  produces  the  minimiom  ripple  at  the  output. 
However,  in  the  case  of  a superconducting  alternator  there  are  some 
potential  problems  associated  with  the  use  of  field  control  for  voltage 
regulation  (see  [10]).  One  area  of  concern  is  the  fact  that  abrupt 
variations  in  field  current  may  cause  this  winding  to  leave  the  super- 
conducting mode.  Another  problem  is  the  long  time  constant  of  the  field 
circuit  which  produces  a very  slow  step  response. 

The  use  of  a phase  controlled  rectifier  bridge  on  the  output  should 
1-rovide  relatively  lower  field  current  variations  and  a faster  step  re- 
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sponse.  It  also  should  be  noted  that  gated  devices  may  be  required 
for  the  rectifier  bridge  in  order  to  provide  short  circuit  protection. 
Thus  these  same  devices  can  be  used  to  provide  voltage  regulation. 


2.2  Steady  State  Model 

The  schematic  diagram  for  the  alternator  with  a controlled  recti- 
fier bridge  is  shown  in  Figure  10,  and  the  output  voltage  and  current 
waveforms  are  shown  in  Figure  11.  The  approximations  used  in  this 
section  are  identical  to  those  used  for  the  analysis  of  the  uncontrolled 
rectifier  bridge. 

2.3  Steady  State  Equations 

Equations  (44.)  - (48.)  of  the  previous  section  provide  five  expres- 
sions which  can  be  solved  numerically  to  find  the  variables  B,  y, 

V and  W.  However,  with  the  thyristor  bridge,  comn.  .tation  does  not  start 

when  V = V,  but  at  some  later  time  when  v > v,  , i.e.,  B is/ controlled 
c b c — b 

externally  to  produce  the  desired  This  means  that  (47.  ) no  longer 

applies  since  it  is  based  on  v^  = Vj^  at  0 = B. 

Therefore,  there  are  only  four  equations  to  work  with, (44.),  (45.), 
(46.)  and  (48.),  and  is  held  constant  slightly  above  the  minimum 
value  requii  -^d  for  maximum  I, . 

Li 

Let,  x'  = 
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6 

y 

V 

w 


(72.  ) 


J 


Figure  11. 


Output  voltage  and  armature  currents  for  the  alternator 
with  controlled  rectifier  bridge. 
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f (x')  = R.H.S.  of 


(44.  ) 
(45.  ) 
(46.  ) 
(48.  ) 


(73.  ) 


A solution  for  x'  can  be  obtained  from  the  standard  Newton-Raphson 
equation, 

F(x')  (x'-  x’)  = f (x')  - f (x') 

—0—0  — — — — o 

where  is  some  initial  starting  point  sufficiently  close  to  }<’  and, 


(74.  ) 


F(x' ) 
— o 


the  Jacobian  matrix  = 


9f ’(x' ) 


(75.  ) 


— — o 


A satisfactory  value  for  ^ is  obtained  by  using  the  previous  diode 
bridge  solution,  jc,  where  is  a variable,  or  by  using  the  approximate 
method  reported  by  Franklin  in  [17,18]. 


2.4  Numerical  Results 


This  numerical  example  is  based  on  the  same  4 pole,  400  Hz.,  10 
MVA/5kV  superconducting  machine  used  in  the  previous  example.  The  same 
full  load  conditions  are  assumed, 


V = 6760  V.D.C. 

IJ 

I = 1420  A.D.C. 

L 


(76.  ) 


The  minimum  value  of  I,  that  will  produce  I = 1420A  corresponds  to  the 

I L 

solution  for  the  uncontrolled  rectifier  bridge  (i.e.,  minimum  8).  This 
current  is  found  from  the  first  section. 


I 


f (min . ) 


250A. 


(77.  ) 
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In  an  actual  system  it  will  be  desirable  to  set  1^.  at  some  value 
above  that  given  by  (77.).  This  will  insure  against  low  voltage  if 
decreases  for  any  reason.  For  this  example,  the  high  value, 
was  arbitrarily  taken  to  be. 


I..  . = 1.1  I.,  . , (78.) 

f(max. ) f(min.  ) 

A plot  of  vs.  for  both  the  controlled  and  uncontrolled  rectifier 
bridge  is  shown  in  Figure  12. 

8,  u,  the  first  three  harmonics  of  v , and  the  rms  value  of  i,  are 

o f 

plotted  vs.  I in  Figures  13  through  18  respectively.  Note  that  the 
L 

corresponding  values  for  the  diode  bridge  case  are  included  for  reference 
purposes . 


i 

f 

As  would  be  expected.  Figure  13  indicates  that  the  thyristor  bridge 
6 must  exceed  the  diode  bridge  6 to  compensate  for  the  higher  I^.  The 

thyristor  bridge  3 also  drops  as  increases  in  order  to  compensate 

t 

for  the  higher  voltage  drop  across  the  armature  windings. 

Figure  14  reveals  an  interesting  characteristic  in  that  the  p for 

the  thyristor  bridge  is  considerably  less  than  the  p for  the  diode 

bridge.  This  is  not  too  surprising  when  Figures  10  and  11  are  considered. 

Referring  to  these  figures,  it  is  observed  that  because  of  the  delayed 

3,  commutation  from  b to  c does  not  start  until  v,  > 0.  Therefore,  more 

, 1 be 

voltage  is  present  to  force  the  commutation  process  than  in  the  diode 
case  where  commutation  begins  at  v =0.  Because  of  this  higher  v , 

DC  DC 

ij^  will  be  driven  to  zero  in  less  time,  thus  producing  a smaller  p for 
the  thyristor  case. 
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Figure  13.  6 vs.  for  both  the  thyristor  bridge 

and  the  diode  bridge. 


3i 


I (amps) 

1j 


Figure  14.  u vs.  I for  both  the  thyristor  bridge 

Li 

and  the  diode  bridge. 
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The  voltage  harmonics  shown  in  Figure  15  to  17  are  obtained  by 
calculating  the  Fourier  coefficients  of  the  output  voltage,  v^,  which 
can  be  obtained  from  (65.)  - (71.).  Again  as  expected,  those  figures 
generally  indicate  higher  v^  harmonics  for  the  thyristor  case  than 
for  the  diode  case.  It  is  also  noted  that  the  harmonics  for  the  thyri- 
stor bridge  tend  to  be  higher  for  the  lower  values  of  I^.  This 
characteristic  is  caused  by  the  higher  6 values  under  light  loading 
conditions . 

The  rms  value  of  i^  can  be  found  by  numerical  integration  of  (27.). 
Figure  18  indicates  that  the  higher  thyristor  8 will  lead  to  higher 
values  of  i^  (rms)  over  the  load  range.  As  was  noted  in  the  discus- 
sion for  the  uncontrolled  rectifier  bridge,  these  i^,  (rms)  calculations 
may  be  hiqher  than  the  actual  values  due  to  a failure  to  account  for  the 
high  frequency  attenuation  of  the  damper  shield. 
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Peak  value  of  harmonic  (volts) 


Figure  16.  12th  harmonic  of  v vs.  I,  for  both  the 

o L 

thyristor  bridge  and  the  diode  bridge. 
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Peak  value  of  harmonic  (volts) 


3.  FAULT  CURRENT  CALCULATIONS 

3.1  Introduction 

This  discussion  is  based  only  on  the  controlled  rectifier’  bridge 

configuration  shown  in  Figure  10.  The  reason  for  this  choice  is  that 

this  circuit  has  the  capability  of  fast  turn  off  in  the  event  of  a 

fault.  Fast  turn  off  can  be  achieved  by  the  system  in  Figure  1 only  if 

some  type  of  series  switch  is  added  to  the  circuit. 

In  addition  to  filtering  the  output  voltage,  the  inductor,  L , in 

o 

Figure  10  must  be  capable  of  limiting  I in  the  event  of  a short  across 

Li 

the  load.  The  length  of  time  that  L^  must  perform  this  fimction  is 

limited  however,  since  the  bridge  can  be  turned  off  on  the  next  cycle 

after  the  fault  is  detected.  It  is  also  common  to  select  L to  limit  the 

o 

C^  charging  current  when  the  system  is  initially  turned  on.  This  can 
also  be  accomplished  by  a further  delay  in  6 however,  so  charging  current 
will  not  be  used  as  a constraint  in  this  analysis. 

To  determine  if  L^  is  of  adequate  size,  it  will  be  necessary  to  cal- 
culate the  transient  load  cu’.’rent,  i.p,  that  occurs  after  the  fault, 
since  the  differential  equations  involved  have  time  varying  coefficients, 
a numerical  solution  will  be  required.  In  this  particular  study  it  is 
assumed  that  the  fault  occurs  at  the  beginning  of  a conduction  period 
and  that  the  bridge  will  not  be  turned  off  until  this  conduction  period, 
the  next  commutation  period,  and  a final  conduction  period  are  complete. 
This  corresponds  to  the  interval,  AE,  shown  in  Figure  19.  The  rationale 
behind  this  choice  is  that  some  time  is  required  for  ij^p  to  exceed 

I,,  V , at  which  time  a current  overload  sensing  circuit  is  enabled  to 
Lv  Tndx  • ) 
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A:  Fault  occurs  at  the  start  of  a conduction  period 
B:  On-coming  thyristors  fire  and  commutation  begins 
C:  Commutation  interval  ends 

D:  I,  , , is  exceeded,  next  trigger  signal  is  blanked 
E:  i^P  decays  to  0 


Figure  19.  Transient  Behavior  for  the  Controlled 
Rectifier  Bridge. 


44. 


blank  all  the  thyristor  gate  signals.  Othei  choice.s  are  certainly  f>os- 
sible,  such  as  assuming  that  the  fault  occurs  at  the  start  of  a commu- 
tation period  and  that  the  thyristor  gate:;  are  blanked  before  the  next 
firing  pulse. 

3.2  Circuit  Model  and  Equations 

Figure  20  represents  the  equivalent  circuit  with  phase  "a"  conducting 

and  a short  across  the  load.  Note  that  the  armature  resistance  R and 

a 

the  parasitic  resistance  of  the  filter  choke,  R , have  been  included 

P 

since  they  aid  in  limiting  the  fault  current.  The  periods  AB,  BC,  etc. 
in  Figure  19  will  correspond  to  the  following  thyristors  being  on: 

Period  Thyristors  Conducting 


AB, 

conduction 

Ql,  Q6 

BC, 

commutation 

Ql,  Q6 

CE, 

conduction 

Ql,  Q2 

The  steady  state  analysis  hac  assumed  that  the  winding  resistances 
are  negligible.  That  practice  will  be  continued  here  for  all  windings 
except  those  on  the  armature;  as  before,  it  implies  that  Xf,  Xd  and  Xq 
are  approximately  constant  over  the  relatively  short  transient  period, 

AE , and  that  the  voltages  across  the  closed  f,  d and  q windings  are  appro- 
ximately zero. 

The  equations  for  the  ce^mmutation  period,  BC,  are  found  to  be, 

di_ 

lAj  ^ = [B]  (79.  ) 

whiere 
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Figure  20.  Equivalent  circuit  for  conduction  and  commu 
tation  periods  while  fault  is  present. 


• md  th»'  |A]  and  [Hi  inati'ices  are  d#*fined  in  Appendix  11. 


The  equations  for  the  conduction  pei i ■!,  AH,  ( Q2  off)  can  be  ex- 


pressed in  a similar  manner, 

'bi'  % = '“ii'i' 

where  sufematrices  of  [A],  [Hi  and  \ These 

submatrices  are  also  defined  in  Appendix  II. 

Since  the  fault  is  assumed  to  occur  at  A in  Figure  19,  (81.)  will 
be  solved  first,  then  (79.).  It  is  unnecessary  to  formulate  equations 
specifically  for  the  second  conduction  period,  CE,  since  this  period  can 
be  analyzed  by  using  (81.). 

As  predicted  earlier,  [A]  and  [B]  have  time  dependent  elements. 

This  implies  t*^  "'.t  i must  be  found  by  some  numerical  integration  tech- 
nique. The  modified  Euler  method  was  chosen  for  this  pai’ticular  study, 
but  other  techniques  could  also  be  employed. 

Starting  with  the  conduction  period,  AB,  (61.'  can  ue  solved  for 
dj_’ 

— it  each  time  increment.  At,  and  the  next  value  of  i'  can  then  be 
dt  - 

found  by  using,  the  .standard  modified  Euler  equations  (see  [27.]).  This 

process  is  simjjlified  however,  if  the  ip[  roxirr.ation  of  constant  Xf,  Xd, 

and  Xq  i.s  again  utilized.  Writing  the  flux  linkage  equations  for  the 

most  general  case  (the  o(jmmut  it  ion  period,  BC),  one  obtains. 


Xf 

Xd 

*1  :o:i 

/3  JOS  (ujt+iT/i)) 

^LF  " 

/J  Mj  sin  (uit) 

f3  M,  sin  (wt) 
d 

Xq 

✓ i '.wtTi/6) 

<1 

/3  M cos  (wt) 
d 

:ont . 


It  7. 


J 


"f  ^fd  ° 

^f(tot) 

L,  0 

i , 

fd  d 

d 

0 0 Lj 

i 

d 

q 

_ _ 

^ 

or  in  vector  notation,  X^,  = xi, ^ + y i,  „ + [C]  i^,  . 

— fdq  — LF  kF  — fdq 


(83.  ) 


—fdq  found  initially  by  substituting  the  steady  state  values  for 

ifdq’  ^LF  ^kF  ~ B+m-h/3  (i.e.,  is  found  from  (25.),  (26.), 

(27.)  and  i,  _ = I,  ,i,  „ = 0).  At  each  time  increment  i,  ^ and  i,  can  be 
LF  L>  kF  LF  kF 

found  from  the  modified  Euler  equations  while  i^,  can  be  found  from  (83.) 

—fdq 

(using  the  constant  value  of  X^,  ). 

—fdq 

Figure  19  indicates  that  the  firing  angle  is  delayed  before  the  fault, 
but  not  afterwards  (point  8).  The  reason  for  this  is  that  once  the  out- 
put is  shorted  the  voltage  regulator  will  call  for  the  minimum  firing 
angle,  and  the  oncoming  thyristors  will  conduct  as  soon  as  possible. 

Since  the  load  current  is  no  longer  constant,  the  steady  state  equations 
that  yield  6(min. ) (i.e.,  the  firing  angle  for  an  uncontrolled  rectifier 
bridge)  do  not  apply,  and  the  firing  angle  (point  B)  must  be  found  by 
determining  the  first  point  at  which  ‘ becomes  positive 


Q2  will  start  to  conduct  and  the  coimutation  of  Q6  will  commence.  V, 
for  the  AB  conduction  state  can  be  readily  found  from  the  be  loop, 

be  a LF  a a dt  f dt  d dt 

_ di 

- /3  M,  cos  (wt)  - 1/3  u)  cos  (cjt)  (M  i + M.  i.) 

d dt  f f(tot)  d d 


be 


+ /3  u)  M , sin  (wt)  i 
d q 


(84.  ) 
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V.  is  tested  at  each  time*  increment  of  the  AB  interval;  at  the  first 
be 

point  where  0 the  computer  program  br'anches  to  tfie  equations  for 

the  BC  commutation  interval.  It  will  also  be  necessary  to  perform  some 
test  to  determine  when  the  commutation  period  ends  at  point  C.  This 
is  done  by  comparing  i^j,  and  ij^^  at  each  time  increment  of  the  BC  inter- 
val; the  commutation  period  ends  at  the  first  point  where  ij^^  ^ 

It  should  be  noted  that  it  may  take  several  machine  cycles  to  com- 
mutate the  thyristor's  after  the  firing  signals  have  been  blanked.  This 
can  be  illustrated  conceptually  by  the  simplified  model  shown  in  Figure 
21  (this  model  is  of  little  quantitative  use  however,  since  it  does  not 
account  for  the  winding  resistances  and  inductances  of  the  machine).  If 
e^^^j  = v^j^  and  the  thyristors  start  to  conduct  at  mt  = tt/3  (approximate), 
i^^j  will  have  the  form. 


/2V 


■(t)  ■ Ig  (1/2  - cos  (a.t)) 

o 


where  I = load  current  at  tut  = tt/3, 
h 


(85.  ) 


At  wt  = ^ (approximate),  conduction  switches  from  phase  b to  phase 

c so  the  source  voltage  becomes  e . , = v . The  load  current  now  has 

s(t)  ac 

the  form. 


L..,  = I,  + — p-  (tt  - cos  (u)t  - 2it/3)) 

( t ) L wL  2 

o 


(86.  ) 


Even  if  the  bridge  is  blanked  duri  g the  v^^  cycle  (to  prevent  the  next 
commutation)  (86.)  never  goes  negative,  meaning  that  the  conducting 
tliyristors  will  not  turn  off. 
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A similar  phenomenon  can  occur  in  the  actual  [physical  circuit, 
except  that  will  eventually  decay  to  zero  due  to  resistive  damp- 

ing. This  may  require  several  cycles  however,  and  more  cycles  will  be 
required  for  larger  values  of  L^.  Thus  if  is  large, several  cycles 
may  be  required  to  complete  turn  off;  however,  if  it  is  too  small  the 
peak  fault  current  will  be  excessive. 

3.3  Numerical  Results 

The  transient  analysis  algorithm  can  be  used  to  plot  post  fault 
current  waveforms  for  various  values  of  L and  R . Two  parametric 

op 

studies  of  i,_  for  variations  in  L and  R are  shown  in  Figures  22  and 
LF  op  ^ 

23  respectively.  Figure  24  shows  a plot  of  i that  requires  three 

Lit 

cycles  for  i^^  to  reach  zero.  Presumably  commutation  would  occur  on  the 

fourth  cycle  where  ij^p  would  attempt  to  go  negative. 

This  algorithm  could  also  be  used  to  plot  during 

a fault  condition.  However,  as  stated  earlier  for  the  steady  state 

calculations,  the  i,,^  values  may  be  too  high  since  the  model  does 
f(tot) 

not  include  the  high  frequency  attenuation  of  the  damper  shield. 
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Figure  22.  Fault  current  vs.  0 for  different  values 

of  L . R = O.OSfi 
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Pigare  23.  Fault  current  vs.  0 for  different  values 

of  R . L =3.0  mH. 
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G (radians) 


figure  24.  Exponential  decay  of  i^^ 
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4.  VARIATION  OF  THE  AETERNATOR  PAPAMETEPS  TO 
DECREASE  OUTPUT  RIPPLE  VOLTAGE 


4.1  Introduction 

As  noted  previously,  there  are  two  basic  methods  for  regulating  the 
dc  output  voltage,  V , of  the  rectified  alternator: 

1.  Use  an  uncontrolled  rectifier  bridge,  and  regulate  V by 

iu 

controlling  the  average  field  current,  I^. 

2.  Hold  I^  constant,  and  regulate  the  voltage  by  means  of  a 
controlled  rectifier  bridge. 

The  first  method  provides  the  minimum  ripple  voltage,  but  it  tends  to  have 
a slow  response  time  due  to  tiie  long  time  constant  of  the  field  winding. 
Therefore  the  analysis  of  this  section  is  based  on  the  controlled  rectifier 
bridge . 

If  the  alternator  is  modeled  by  an  ideal  ac  voltage  in  series  with 
an  inductor,  it  is  well  known  that  an  increase  in  ■t'his  inductance  will 
increase  the  commutation  angle,  h,  shown  in  Eigui’e  11.  For  constant  V 

L 

and  I,  this  effect  can  also  lead  to  a reduction  in  output  ripple  voltage, 

Lj 

as  illustrated  in  Figure  25.  Comparing  parts  (a.)  and  (b.)  of  the  figure 
it  is  seen  that  the  same  average  output  voltage,  V , is  achieved  by  dif- 
ferent  combinations  of  the  angle,  d and  p.  However,  the  deviation  of 

V is  less  in  part  (b.).  This  indicates  than  an  iru.re.ije  in  ij 
o 

requires  the  bridge  firing  angle  to  decrease,  resulting  in  a more  level 

V . Therefore,  the  R„,  p„  combination  in  (b.  ) produces  a lower  ripple 

o 2 2 

voltage  than  the  combination  in  (a.).  This  imi>lies  tliat  at  least 

over  a limited  range  of  p values  it  is  possible  to  deci'oase  the  rip;  le 
voltage  by  increasing  p.  Tiius  for  a fixed  load,  it  is  possible  to  decrease 
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average  output  voltage 


(source  inductance  for  (b 
higher  than  that  for  (a.) 


Figure  25.  Effect  of  p upon  output  ripple  voltage  for 


. ) is 

) 


fixed  V,  and  I . 


the  ripple  by  increasing  the  source  inductance  since  this  causes  an 
increase  in  ij. 

The  model  used  in  this  study  was  considerably  more  compl icateii 
than  the  one  just  described,  but  it  was  found  that  for  a constant  load 
the  ripple  could  be  decreased  if  the  armatui’e  inductance,  L , was 
increased  beyond  its  specified  value  of  0.3  mH.  This  implies  that  a 
ligtiter  weight  filter'  could  be  used  if  L were  increased.  A word  of 

3 

caution  is  in  order  however’,  since  these  gains  m-iy  b--  iffi-.'t  t.y  an  in- 
crease in  alternator  weight  (due  to  the  1 irget’  L ).  bltimatoly,  it 

3 

would  be  desir.tble  t<’j  develoj'  a piw  edura-  tii.it  would  f.inin.i^e  th»-  ■ om- 
bined  weight  of  the  alternator  and  liltei  . This  w...uM  re.pairt  a detailed 
weight  analysis  of  the  su[>erconducting  alternator  however,  ml  :u<'h  an 
effort  would  be  beyond  the  scope  of  this  pre.ierit  stu.lv. 


4.2  Effect  of  on  the  Output  Voltage  Harmonics 

As  discussed  in  the  previous  section,  it  is  possible  to  reduce  the 
full  load  ripple  voltage  by  changing  the  output  impedance  of  the  alter- 
nator. The  model  used  in  this  study  cannot  be  reduced  to  a single  ac 
source  in  series  with  such  an  impedance,  but  a similar  effect  will  occur 
if  the  armature  self  inductance,  L , is  varied.  Changes  in  L will,  of 

3 3 

course,  change  the  armature  mutual  inductances.  To  account  for  this, 
it  is  assumed  that  the  coefficien""  of  coupling  between  and  each  of  the 
other  windings,  remains  constant  wtiile  L is  varied,  i.e., 

3 

”f 

k ^ = = constant  ( 87 . ) 

df  rr~r~ 
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constant 


(88.  ) 


constant 


(89. ) 


(L^  and  L,  also  remain  constant.) 
r d 


4.3  Numerical  Results 

Since  )<  , and  are  assumed  constant,  (87.)  indicates  that  an 
af  f 

increase  in  L will  also  increase  M^.-  Thus  an  increase  in  L implies 
a f a ^ 

that  the  same  magnetic  flux  linkage  from  field  to  armature,  can 

be  achieved  with  a lower  (since  a thyristor  bridge  is  used  for 
voltage  regulation  it  is  assumed  that  will  be  held  constant  at  110% 
of  the  minimum  allowable  value  for  a given  L , as  discussed  in  the 

cl 

section  on  controlled  rectifiers).  Figure  26  indicates  the  decrease 

in  the  required  as  is  increased  for  = 1420  A.  dc.  This 

decrease  in  might  allow  the  use  of  smaller  superconducting  wire  for 

the  rotor  winding,  thus  decreasing  the  rotor  size.  An  alternate  appr  ach 

would  be  to  hold  constant  and  allow  L^-  to  decrease  as  L increased; 

f f a 

thus  the  field  winding  would  have  fewer  turns  (both  effects  appear  small). 

Figure  27  indicates  that  the  6th  harmonic  of  v reaches  a minimum 

o 

at  L = 0.72  mil.  This  leads  to  a decrease  in  the  size  of  the  output 
filter,  since  less  attenuation  is  required. 

Break  frequency  = f = — — (90.) 


Figure  27  also  indicates  that  the  12th  and  18th  harmonics  generally 


coi  tinue  to  increase  with  L , but  their  effect  is  less  important  since 

3 
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filter  attenuation  is  much  greater  at  these  frequencies.  F’igure  28 

shows  that  the  6th  harmonic  continues  to  decrease  as  increases, 

indicating  that  the  optimum  at  ^0%  of  full  load  lies  somewhere 

above  0.9  mH.  Thus  the  optimum  L is  different  for  different  values 

^ a 
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ricure  27.  6th,  12th  and  18th  harmonics  of  v vs. 
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for  the  controlled  rectifier  bridge. 

I,  = 1420  A.d.c. 

L 


61 


MIHIMTZATION  OF  L C FILTEF  WEIGHT 
o o 

5.1  Introduction 

The  previous  sections  have  considered  the  following  topics: 

1.  Steady  state  behavior  with  an  uncontrolled  rectifier  bridge. 

2.  Steady  state  behavior  with  a controlled  rectifier  bridge. 

3.  Transient  cur-rcnts  that  occur  when  a short  circuit  is  placed  across 
the  output. 

4.  Reducing  the  full  load  output  ripple  voltage  by  increasing  L . 

d 

The  results  of  these  studies  can  now  be  used  in  designing  an  output 

filter  for  minimum  weight.  This  analysis  assumes  the  use  of  a controlled 
rectifier  bridge  for  voltage  regulation.  The  maximum  allowable  ripple 
voltage  will  be  based  on  the  size  of  the  sixth  harmonic  that  is  present 
at  full  load.  It  should  be  noted  that  this  harmonic  will  actually  be 
greater  at  minimiim  load  since  the  firing  angle  of  the  thyristors  will  be 
greater.  This  study  assumes  that  the  load  will  b-.  fairly  constant  how- 
ever, and  that  the  presence  of  ripple  voltage  will  be  more  important  at 
full  load  than  at  lighter  loads.  Hence  the  filter  optimization  is  based 
on  full  load  conditions. 

5.2  Calculation  of  L and  C for  Minimum  Total  Filter  Weight 

o o 

In  the  weight  minimization  algorithm,  and  are  calculated  to 
provide  a given  amourit  of  ripple  ••ttenuation  at  full  load.  Tliis  calcula- 
tion is  based  only  on  the  sixth  harmonic  and  ignores  the  harmonic  attenua- 
tion provided  by  the  load  in  conjunction  with  L^.  Therefore, 
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1-L  C oj 
006 


(91  . ) 


where 

and 


specified  magnitude  of  the  6th  harmonic  attenuation 


li)  . = 15079.64  radians/sec. 

b 


L C > 
00  — 


A!i 

2 

k ui 
16 


(92. ) 


Due  to  the  high  value  of  the  magnetic  field  and  the  low  weight 
requirement,  it  is  assumed  that  will  be  an  air  core  reactor. 

Aluminum  was  chosen  for  the  conductor  due  to  its  low  weight/conductance 
ratio.  The  physical  configuration  of  the  inductor  is  shown  in  Figure 
29.  This  particular  design  is  chosen  to  produce  the  minimum  loss  for 
a given  amount  of  material  (see  [25,26]).  The  inductance  is. 


L = (24. 5 X 10‘^)  a H.  (93.  ) 

o 

cind  the  weight  of  is  given  by 

L Wt.  = 3 TI  f D a^  lbs.  (94.  ) 

o w 

I where 

N = number  of  turns 

f ■ 

a = thickness  of  the  coil  (m.  ) 

' D = density  of  A1  (5837.8  lb‘/m.^) 

w 

1 f = filling  factor  of  the  conductors  (assumc-d  to 

be  0.7) 

The  filling  factor  can  be  expressed. 
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Figure  29.  L dimensions. 


f 


(95.  ) 


NF 

w 

2 
a 

2 

where  F = cross  sectional  area  of  one  windineCm.  ) 
w “ 

The  common  method  of  specifying  capacitor  weight  is  in  terms  of 

jouls/lb.  Therefore  the  total  weight  of  the  L C filter  can  be 

o o 

expressed. 


T = C weight  + b weight 
wt  o o 


C V 

- + 3 Tt  f D a^  lbs. 

2D  w 


where  D = energy  density  of  C ( joules/lb . ) . 
c o 

Therefore  substituting  (92.),  (93.)  and  (95.)  into  (96.) 


bt  = S'*'  ' S*' 


where  = I 


VSl.lxlO 


-7 


F V \ / k +1 
-2LiiW  _i 


Vd  k w 

' c 1 c 


k^  = 3 Ti  f D 

3 w 


To  find  the  minimum  value  of  T set 

wt 


dT  ^ -5k„  _ 

^ . 3 k = 0 

da  6 3 


5k^  \l/8 

I 


Once  a is  determined,  N can  be  found  from  (95.)  , L from  (93.) 

o 

and  C from  (92.  ) . 
o 


(96.  ) 
(97.  ) 


(98.  ) 
(99.  ) 
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5.3  L C Design  Algorithm 
o o 

The  flow  chart  for  the  L C filter  design  program  is  shown  in 

o o 

Figure  30.  Tlie  following  discussion  refer’s  to  the  seven  blocks  indi- 
cated in  the  figure. 

1.  Read  input  data.  This  includes  the  following  information: 

etiergy  density  (D^),  current  density,  maximum  6th  harmonic  ripple 
voltage  at  full  load,  and  maximum  short  circuit  current,  I,,  ,. 

This  particular  progi’am  assumes  that  the  following  quantities  are 
constant : 

V,  = 6760  V.dc 

L) 

I,  = 1420  A.dc 

Lj 

= 15079.6  rad. /sec.  (line  frequency  = 400  Hz.) 

6 

V,  , I,  and  01.  could  be  varied  if  desired,  by  making  a few  minor  changes 
L L 6 

in  the  program. 

2.  Set  L = 0.3  mil,  the  normal  design  value  specified  in  [13,14]. 

3.  Assuming  that  a controlled  rectifier  bridge  is  used,  the  minimum 
weight  and  that  will  meet  the  6th  hannonic  ripple  specification 
are  calculated. 

4.  A transient  analysis  subroutine  to  called  to  determine  if  the  peak 

short  circuit  current  will  exceed  the  specified  value  of  I,  , Hie 

details  of  this  analysis  are  given  in  section  3. 

5.  If  I,  / s is  exceeded,  L Ls  increased  by  10%,  and  step  3 is  repeated 

L(max.)  o i ■>  t 

(C  is  simultaneously  decreased  to  maintain  a constant  L C product]  If 
o o o ' 

I,/  \ is  not  exceeded,  the  L and  C design  data  is  printed. 

L(max)  o o ‘ 

6.  The  program  calculates  and  C^,  first  for  the  normal  and  then  for 
the  optimum  values  of  I.  . If  the  calculation  for  the  optimirm  L,  has 
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been  completed  the  program  ends.  If  not,  the  program  branches  to  step  7. 
7.  Set  L = 0.72  mil.,  the  optimum  value  toi'  minimurii  r-ipple  at  full 
load  calculated  in  Section  4.  Steps  3 through  6 are  then  repeated. 

5.4  Numerical  Res<ilts 

A sample  of  the  computer  results  for  the  optimization  program  are 
shown  in  the  following  example.  Note  that  use  of  the  optimum  L 
decreases  the  total  filter  weight  by  approximately  22%. 

The  total  filter  weight  will  obviously  decrease  if  a higher  energy 
density  is  used  for  and/or  a higher  current  density  is  used 

for  L^.  Plots  of  filter  weight  vs.  energy  density  and  current  density 
are  shown  in  Figures  31  and  32  respectively.  The  filter  weight  will  also 
be  affected  if  the  allowable  I,,  \ is  changed.  A plot  of  this  is 

shown  in  Figure  33. 
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VriT^  THE  F0LL0V'It.JG  PAf^^’IETERS  F-if.  THE  METER 


ALL  inPUTS  HAVE  FDRMAT  = F7 . 2 UHLESE  OTHERWISE  SPECIFIED 
CA"’.  "'•H’^GY  DE.JSITY  C H JT'LE.c/lb.  ) = 

bn.  0 


Cl'PP.l.’T  PE.3SIT';  PEP  E "IPE  (CIP  .,IL/A.IPI  - 

'3.  G 


ha::.  p;:s  -'aeve  ep  =.th  ;ah?i.  ji  vo  (V..lts)  = 
■D.  C 


,»LLJ'..EP,LE  PEAK  FAIG.T  CURREHT  CA.IPS)  = 

Dsno . n 


THL  FCLLC-’IHG  VALI:ES  APE  BASED  Gi:  rJJFIE'U.  LA  = C.30v3E-Q3  H. 

IF=  0.27SE+02  EETA=  0.212E+01  :1U=  0.  307E+00 

IL=  0.1A2E+0A  7L=  0.67  3E+OA  Vl^  O-IHOE+OA 

0PT.7,\LUE5  EEFPRE  FALLT  TEST  ARE  LO=  0.A22E-02  H..  C0=  0.3C6E-0A  ]•  D . 

ES=  0.5A7E-01 

MAX  L0AD  CURREHT  FP0M  FAULT  = 0.  250E+0A 

FAULT  CURRSOT  T00  LARGE,  LO  IH  CREASED 
MAiX  L0AD  CURRENT  FROM  FAULT  = 0.  250E+0A 

MAX  L0AD  CURRENT  FR0M  FAULT  = 0.2A5E+0A 

L3=  0.511E-02.  CO=  n.r.23E-0A  ILF--  0.  2A  f E'- ''A  REG=  '^.C-ISE-Pl 

L"T=  r.'‘42L+r,2  C"T=  C.2Ai3E+02  TCTmL  "T  - P.  7 ' M-*-QZ 

.....  TE:e;J  = l 3-'l  L radius  = n.  223L+0PM  L ..EiGTM  = C . 1 1 1 S + ' 

T:HE  FALL.iWItJ  i VALUES  APE  BASED  ONaPTIMUM  LA  = p.  7^0E-  33 
IF=  Q.233L+03  PETA=  0.223E>01  PU=  0.6uA_+00 

IL=  0.1/12E+CA  7L=  0.376H-CA  71=  0.620E+n3 

CPT.7ALUES  •:EF;PL  FALT.P'  test  ARE  L0=  0.2';2L-C.P  E.,  C0=  0..',3'>E-~-  ' 

ES=  0.A23L-0  1 

MAX  LOAD  CURREiJT  FP.J;;  FAILT  = 3.  250E+0A 
;,.UL7  CUP"E.IT  TG.;  L-\P3E,  LG  I3CPLASLL 

;:.v:  LcAL  cl;;pe.jt  fe....  fault  = 3.  250E+o^. 

MAH  LOAD  CurPE-JT  FROM  FAIET  = C.  250E+04 

EG:  -.',a,1L-03  C3=  0.A1D.--3A  ILF=  3.  S E..-'- A P.:,S=  1 

L-’T=  P.4P5E+0'’  C”T=  0.133L+02  TU  T.'tL  V'T_^  0.613E+32 

. TU':;p  = 11/  L FAD.1H3  = o.  206E+on;  E l_.ht:;  : o.l03_^•  3.; 

•..'PITE  "3"  7.,  E.JD,  .M'.  "1"  F.FP  AMJTHER  RUi; 
f 0;iMAT=  I 2 

0 
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CAPACITOR  DENSITY  (Joules/lb.) 


f'ii;ure  31.  Filter'  weight  vs.  energy  density.  Current 
den.sity  = 100  cir.  mile/amp,  dll  other 
variables  arc  the  same  as  in  the  example  run. 


6.  SENSITIVITY  ANALYSIS 


6.1  Introduction 

Since  the  values  of  the  alternator  inductances,  L , L^,  L.,  M , M_, 

^ i Q d I 

M,  and  are  subject  to  numerical  error,  it  is  of  interest  to  see 

d fd 

how  errors  in  these  parameters  will  affect  the  calculations  for  I^,  6, 
y,  V and  W.  As  noted  in  Sections  1 and  2,  the  calculations  for  the  uncon- 
trolled and  controlled  rectifier  bridges  are  quite  similar.  This  implies 
that  the  effect  of  a given  parameter  error  should  be  about  the  same  for 
both  types  of  systems.  Therefore  it  was  decided  to  limit  the  sensitivity 
analysis  to  the  uncontrolled  rectifier  case.  For  convenience  we  define 


(100. ) 


the  following. 


V 

W 


Z = 


(101.  ) 


74. 


(44.  ) 


f^(x,y)  = R.H.S.  of 


Theoretically,  this  analysis  could  be  performed  by  either  of  two 
methods : 

3x 

1.  Use  (102.)  to  find  ^ and  solve  for  Ax  for  a given  A^, 

i.e.  , Ax  = Ay.  This  will  be  referred  to  as  the  differential 
- 

method. 

2.  Simply  replace  by  + A^  and  use  the  Newton  Raphson  method 

to  find  the  resulting  x + Aj<.  This  will  be  referred  to  as  the  deliberate 
error  method. 

The  differential  method  is  certainly  the  more  elegant  of  the  two, 
so  this  was  investigated  first.  Unfortunately  thus  approach  depends  on 
solving  sets  of  simultaneous  equations  that  have  ill-conditioned  coeffi- 
cient matrices.  Two  algorithms  were  used  for  solving  these  equations, 
but  both  failed  due  to  excessive  round-off  errors.  Therefore  it  was 
necessary  to  resort  to  the  deliberate  error  method.  This  second  approach 
worked  satisfactorily  even  though  it  is  rather  inefficient  in  terms  of 
computation  time.  Both  methods  will  be  discussed  for  completeness,  even 
though  the  first  did  not  produce  satisfactory  results. 

6.2  Differential  Method 

One  usually  does  not  bother  to  describe  methods  that  do  not  work,  but 
this  analysis  is  interesting  from  a conceptual  standpoint,  so  it  is 

7b. 


(45.) 
(46.) 
(47.  ) 
(48.  ) 


(102. ) 


included  for  that  reason.  It  is  also  possible  that  the  problems  with 
this  method  may  eventually  be  solved,  even  though  it  was  unsuccessful 
in  this  present  research. 

Taking  the  partial  derivative  of  (102.)  produce's  an  equation  of 
the  form, 

— ^ = [C]  ^ t r = 0 (103.) 

3 V . 3 V . — — 

: 1 

where  [C]  is  a (5x5)  coefficient  matrix,  and  r is  a (5x1)  vector. 

-1 

^ = - [C]  r (104.) 

9yi 

3}^ 

which  is  the  ith  column  of  , a (5x7)  matrix.  Therefore  it  is  con- 

^iL 

ceptually  possible  to  use  (104.)  for  all  seven  elements  of  to  find 
3x 

. Ax  for  a given  Ay  is  then, 

3y_  - 

3)< 

Ax  = — Ay^  (10  5.) 

Unfortunately,  the  [C]  matrices  indicated  by  (103.)  are  very  ill 

conditioned  in  this  application,  and  this  prevented  finding  a solution 
3x 

for  — . Two  methods  of  solution  were  attempted,  the  first  being  the  DGELG 
3^ 

double  precision  subroutine  from  the  IBM  Scientific  Subroutine  Package 
and  tlie  second  being  a Shipley-Coleman  inversion  algorithm  to  find 
[Cl  Both  of  these  programsuse  pivoting  for  size,  but  they  were  still 
incapable  of  finding  the  correct  solution.  Therefore  this  approach 
was  abrandoned  in  favor  of  the  deliberate  error  method. 


76. 


6.3  Deliberate  Error  Method 


In  this  method  a given  error,  is  added  to  and  the  result- 

ing ^ + A>{  is  calculated  by  the  equations  described  in  section  1.  The 
terms  of  ^ are  not  independent  however,  since  mutual  inductance  terms 
are  present,  and  this  must  be  accounted  for  in  the  analysis.  There- 
fore, the  approach  used  in  this  particular  study  was  to  assume  that  the 

following  terms  can  be  varied  independently  of  one  another:  L , L_,  L,, 

r 

k , k Cl  k ,,  and  k.,,  where  the  last  four  terms  are  the  coefficients 
aa  af  ad  fd 

of  coupling,  i.e.. 


aa 


af 


/llT 

a f 


ad 


/lTT 

a d 


fd 


fd 

•v; 


(106. ) 


The  independent  and  dependent  parameters  are  listed  as  follows: 


Independent 


Dependent 


aa 


af 


ad 


M^,  M^, 


”f’  "fd 


"d’  "fd 


fd 


fd 


for  example,  a 10%  increase  in  L implies  (new  value  = 1.1  L ), 

a 
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M.  = 1.1  k_  L_,  = k,,  ^1-1  " *^ad  ’ 


aa  a’  f af 


ad 


whereas  a 10%  increase  in  k implies  (new  value  = 1.1  k ,), 

ai  al 


M . = 1.1  k ^ /lTT 

af  af  a f 


The  effect  of  these  errors  will  be  described  in  the  next  section  on 
numerical  results. 


6.4  Numerical  Results 

The  following  paragraphs  discuss  the  effects  of  varying  each  of 
the  machine  inductances,  i.e. , the  effect  of  a deliberate  error.  Note 
that  since  the  algorithm  used  in  Section  1 depends  upon  the  approxima- 
tion given  by  (41.),  it  is  necessary  to  restrict  the  parameter  varia- 
tions to  the  range  where  (41.)  is  valid.  It  is  assumed  that  (41  . ) is 
satisfied  as  long  as  the  following  condition  is  met; 


M - M < 0.1  M 

o oo  oo 


AL  : Results  are  shown  in  Figures  34  and  35.  These  figures  indicate 

that  all  of  the  ^ variables  are  quite  sensitive  with  respect  to  AL  . 

AL^:  Results  are  shown  in  Figures  36  and  37.  It  is  noted  that 

tJ,  p,  V and  W do  not  vary  with  respect  to  L^.  Tlie  reason  for  this  is 
that  the  algorithm  automatically  adjusts  I^  to  comieiisate  for  any 
'hanges,  so  that  the  Mjlj-  flux  linkages  remain  constant  (note  that 
is  dependent  on  L^.)  Compare  with  the  results  shown  in  Figures  42  and 


43. 


o 

0.60 


n 

0.80 


“I 1 

1.00  1.20 

(H.) 


“I 

1 .kO 


Figure  36. 


B and  p variation  vs. 


S- 


n 

1 .60 
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AL^:  Results  are  shown  in  Figures  38  and  39.  These  figures 

indicate  that  the  calculations  are  completely  insensitive  to  varia- 
tions. The  reason  for  this  stems  from  tlie  previously  mentioned  approxi- 
mation , 


M M 
O oo 


(107. ) 


which  is  given  by  (41.)  in  Section  1-  Once  this  approximation  is  made, 

is  replaced  by  in  all  subsequent  calculations.  The  value  of 

M is , 
oo 


oo 


(108. ) 


The  result  of  this  is  that  does  not  actually  appear  in  (44.) -(51.), 
which  are  the  equations  used  to  find  x. 


Ak  (AM  ):  Results  are  shown  in  Figures  40  and  41.  The  figures 

indicate  that  p is  quite  sensitive  to  AM  variations,  while  I,  and  8 
' a f 

are  less  sensitive.  W and  V also  vary  considerably  with  respect  ;o  M^. 


Ak  (AM  ):  Results  are  shown  in  Figures  42  and  43.  For  an  explana- 

3 1 r 

tion  of  these  results,  refer  to  the  discussion  for  AL^. 


Ak^^  (AM^):  Results  are  shown  in  Figures  44  and  45. 


Ak^^  (AM^^):  Results  are  shown  in  Figures  46  and  47.  These  figures 

indicate  that  the  calculations  are  'nsensitive  to  M^,  variations.  The 

fd 

reason  for  this  is  much  the  same  as  for  L^,  i.e.,  once  tlie  approximation 
(107.)  is  made,  M^.^  does  not  appear  in  any  of  tfie  subsequent  eqvjations. 
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L,  (uH.) 
d 

Figure  38.  3 and  m variation  vs.  L 
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Figure  45.  I-,  W and  V variation  vs.  M,. 
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To  suirmarize,  it  can  be  seen  that  g,  yi,  V and  W as  a group 

are  most  sensitive  to  errors  in  L , M and  M,.  is  also  quite  sen- 

a a d f ^ 

sitive  to  errors  in  and  M^,  but  since  the  algorithm  adjusts  to 
maintain  a constant  product  (mutual  flux  linkages)  errors  in 

and  have  virtually  no  effect  on  6,  u,  V and  W. 
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7.  VOLTAGE  REGULATOR  AJ4D  CURRENT  OVERI.OAD  PROTECTION  CIRCUITS 

7 . 1 Introduct ion 

The  experimental  portion  of  this  study  consisted  of  designing  and 
building  a phase  controlled  voltage  regulator  for  eventual  testing  with 
an  alternator.  This  circuit  was  designed  and  tested  early  in  the 
project  when  it  was  thought  that  the  alternator  could  be  modelled  by  a 
voltage  source  in  series  with  single  inductance.  Testing  the  regulator 
with  this  typo  of  a source  would  have  been  a fairly  simple  matter,  but 
such  a test  now  appears  to  have  limited  value  since  a more  detailed  ma- 
chine model  w.is  employed.  Therefore  it  was  decided  to  concentrate  more 
effort  on  the  analytical  study  and  postpone  this  testing  until  a con- 
ventional alternator  with  known  inductances  could  be  obtained.  Schematics 
of  the  complete  design  are  shown  in  I’igures  48  and  49.  Operation  of  the 
voltage  regulator  circuit  in  Figure  48  is  described  in  [23],  and  the 
operation  of  the  current  overload  circuit  in  Figure  49  is  described  in  [22]. 
The  parts  list  is  shown  in  Table  I. 

7.2  Experimental  Results 

As  stated  above,  the  experimental  results  were  limited  to  building 
and  testing  a phase  controlled  voltage  regulator  circuit  with  a current 
overload.  This  circuit  has  the  following  characteristics: 

1.  The  output  voltage  can  be  varied  continuously  from  0 to 
290  V.d.c.  with  a 100  ohm  load. 

2.  The  overload  circuit  turns  the  regulator  off  at  a load  current 
ot  approximately  3.5  A.d.c. 
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TABLE  I.  PARTS  LIST  FOR  CIRCUITS  IN  FIGURES  48  AND  49 


R15  lOKO 

R37  4700 

R7  3 20Kii 

R25  250KO  trimpot 

R38  4700 

R74-P79  300  - 1 watt 

R27  4.7KiI 

R39  4700 

Cl  O.lvif 

R26  3.76Kf2 

R40  4.7KO 

C2  O.OSpf 

R1  IIKSI  - 2 watt 

R41  4.7KO 

C3  lOOuf 

R28  33'in 

R42  4.7K0 

C4  O.lpf 

R29  5.6K« 

R43  4.7KO 

C5  O.OSuf 

R30  250KO  trimpot 

R44  4.7KO 

C6  lOOuf 

R16  1.5Kf2 

R45  4.7KO 

C7  O.luf 

RIO 

R46  100 

C8  lOOpf 

Rll  10K« 

R47  100 

C9  O.OBpf 

R12  250KJ2  trimpot 

R48  100 

CIO  O.OOlpf 

R13  4.7KS2 

R49  100 

Cll  O.lpf 

R14  3.76Kn 

R50  100 

C12  2pf 

R2  11KS7  - 2 watt 

R51  100 

C13  O.lpf 

R9  330n 

R52  3.3KO 

C14  4pf 

RJ7  5.6KSI 

R53  20KO  trimpot 

CIS  O.luf 

R18  250Kfi  trimpot 

R54  2000 

C16  O.luf 

R19  l.SKn 

R55  20KO 

C17  O.luf 

R20  1.5KO 

R56  0.040 

C18  O.luf 

R21  lOKlJ 

R57  lOKO 

T1-T6  Sprague  11Z12 

R22  250K0  trimpot 

R58  20KO  trimpot 

Q1-Q3  2N3906 

R23  4.7KS2 

R59  2.0KO 

Q4-Q9  2N3904 

R24  3.76KfJ 

R60  6300 

Q10-Q15  2N2222 

R3  llKi;  - 2 watt 

R61  l.lKSi 

Q16,Q17,Q18,Q20,Q21,Q22,Q24  2N2222 

R4  330fi 

R62  lOKO 

Q19  2N2646 

R6  5.6HIi 

R63  lOKO 

Q23  2N2907 

R5  250KL  trimpot 

R64  2.2KO 

CR1-CR12  1N914 

R7  lOKii 

R65  10K52  trimpot 

CR13-CR18  1N4001 

R8  lOKfi  trimpot 

R66  15KO 

K1  Reed  Relay  SPST  N/0 

R31  l.SKft 

R67  330K0 

A1-A3  Telefunken  UAA145 

R32  1.5KS2 

R68  330KO 

A4-A5  Fairchild  9601 

R33  l.bidi 

R69  10W2  trimpot 

A6  Signetics  7493 

R34  470n 

R70  15Ki2 

A7  National  Semiconductor  N7408 

R35  470fi 

R71  lOKSi 

SCRl-SCRG  2N1849 

P36  4700 

R72  47KI2 

98. 


3.  No  misfiring  problems  were  observed  once  the  final  design  was 
complete. 

Certain  wavefoi’ms  of  interest  are  shown  in  Figures  50  through  53. 
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8.  CONCLUSIONS 


This  study  indicates  that  it  is  possible  to  utilize  to  help 

perform  some  of  the  functions  normally  assigned  to  the  L C output 

o o 

filter.  This  implies  that  a smaller  filter  can  be  used,  thus  decreasing 
the  weight  of  the  and  C^  components.  For  the  10  MVA/SkV  example 
alternator  with  a controlled  rectifier  bridge  it  was  shown  that  an  in- 
crease in  L from  0.3  mH.  to  0.72  mH  decreases  the  filter  weight  by 
about  17  lbs.,  a 22%  reduction.  This  example  also  indicated  0.72  mH.  to 
be  an  optimum  value,  i.e.,  filter  weight  increased  for  L > 0.72  mH. 

3 

Naturally,  this  savings  may  be  offset  by  an  increase  in  alternator  weight 

due  to  the  larger  L . Therefore  any  final  weight  optimization  study 
3 

should  consider  the  alternator  and  filter  as  a combined  system. 

In  the  course  of  developing  the  filter  weight  minimization  program 
it  was  necessary  to  derive  both  steady  state  and  transient  models  for  the 
alternator  and  rectifier  bridge.  Because  of  the  large  amounts  of  infor- 
mation provided  by  these  models,  it  appears  they  may  be  useful  for 
simulating  the  system  during  the  design  stage.  Once  experimental  data 
becomes  available  for  comparison,  these  models  may  be  refined  as  neces- 
sary in  order  to  accurately  predict  the  various  winding  currents,  com- 
mutation angles,  etc.  It  is  stressed  that  this  experimental  verification 
is  necessary,  and  plans  have  been  made  to  proceed  with  this  for  a system 
with  a conventional  alternator. 
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9.  RECOMMENDATIONS 


This  study  indicates  that  if  is  increased  up  to  a certain 
optimum  point,  it  is  possible  to  significantly  reduce  the  size  of  the 
output  filter.  Information  of  this  type  should  be  brought  to  the 
attention  of  machine  designers,  but  it  may  or  may  not  influence  the 
design  of  future  alternators  due  to  the  many  other  factors  which  govern 
the  size  of  L^.  Ultimately  the  alternator  and  filter  should  be  con- 
sidered together  in  future  weight  minimization  studies. 

Perhaps  the  most  pressing  need  at  this  point  is  to  obtain  some  ex- 
perimental data  to  compare  with  the  predicted  results.  Eventually  this 
must  be  done  using  a superconducting  alternator;  however,  it  is  unlikely 
that  such  a machine  will  be  available  for  this  purpose  for  quite  some 
time.  In  the  interim,  it  is  proposed  that  tests  should  be  conducted  on 
a conventional  alternator-rectifier  system  in  order  to  evaluate  the 
models . 
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APPF;NDIX  I:  GLOSSARY  OP  TERMS 


A,  B,  C = constants  defined  by  (*49.),  (50.)  and  (51.) 
a = thickness  of  L 

o 

= output  filter  capacitor 
= energy  density  of 
= density  of  aliiminum 
f = fill-in  factor  for  L 

o 

fj^  = break  frequency  of  output  filter 

F = cross  sectional  area  of  one  winding  of  L 
w o 

i , i,  , i = line  cui'rents 
a b e 

i^,  i^  = currents  in  the  equivalent  direct  and  quadrature  windings 
used  to  represent  the  damper  shield 

= average  field  current 

i^  = time  varying  component  of  field  current 

i,  = commutation  current 
k 

i^j,  = load  current  with  short  circuit  across  load 

i,  _ = coiranutation  current  with  short  circuit  across  load 

kP 

i^  , i.  , i = field  and  damper  currents  at  0 = B+u-tt/3 
fo  do  qo  * 

K , K,,  K.  = constants  defined  by  (15.)  and  (18.) 
q f d 

k = coefficient  of  coupling  between  armature  phase  windings 

k , = coefficient  of  coupling  between  armature  and  equivalent 

<3  Q 

damper  windings 

k , = coefficient  of  coupling  between  field  and  armature 

af 

kj.^  = coefficient  of  coupling  between  field  and  equivalent  direct 
axis  damper  windings 

kj^  = specified  fiariwuiic  attenuation  factor  of  output  filter 
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k^,  kg  = constants  defined  just  below  (97.) 

I = load  current 
L 

L = self  inductance  of  each  armature  winding 

= self  inductance  of  the  direct  and  quadrature  axis  windings 
= self  inductance  of  the  field  winding 
= output  filter  inductor 

M = magnitude  of  mutual  inductance  between  armature  windings 

3 

= magnitude  of  mutual  inductance  between  damper  and  armature 
windings 

= magnitude  mutual  inductance  between  field  and  armature  windings 

Mfd  = mutual  inductance  between  field  and  direct  axis  damper 
windings 

M , M = constants  defined  by  (33.)  and  (34.) 
o oo 

N = number  of  turns  for  L 

o 

R = resistance  of  one  armature  winding 

3 

R = resistance  of  inductor,  L 
P o 

v^,  V.,  V = voltages  across  rotor  windings 
t d q 

v^,  = phase  to  phase  armature  voltages 

v^  = instantaneous  rectifier  output  voltage 

V = average  output  voltage 

iJ 

V = variable  defined  by  (21.) 

W = variable  defined  by  (20.) 

B = angle  at  which  commutation  starts 
= constant  defined  by  (36.) 

Aj  = constants  defined  by  (36.) 
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^ J ^ J ^ 

a b c f d q 


flux  linkages 


9 = time  angle 

U = commutation  angle 

ui  = electrical  angular  velocity 
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APPENDIX  II:  TRANSIENT  MATRICES 


-/3  M u)  cos  (wt)  -/a  M w cos  (oot)  /s  M cj  sin  (lot)  -2R 


APPENDIX  III;  MAIN  PROGRAMS 


The  following  programs  are  listed  in  alphabetical  order.  All  sub- 
routines except  GELG  and  ARCSIN  are  listed  in  APPENDIX  IV.  It  should  be 
noted  that  the  notation  in  the  programs  occasionally  differs  from  that 
in  the  text: 

Text  Program 


ab 


1.  CONT : Finds  the  solution  for  the  controlled  rectifier  bridge  case. 

2.  MACH2:  Finds  the  value  of  which  produces  the  minimum  value  of 

the  6th  harmonic  of  v . 

o 

3.  MAST : Finds  the  minimum  filter  weight  for  a given  set  of  specifi- 

cations 

4.  PLTDAT : General  purpose  program  that  includes  /arious  simulations  for 

both  the  controlled  and  uncontrolled  rectifier  bridge. 

5.  SENSI3 : Sensitivity  analysis  program. 

6.  TABLE:  Determines  the  harmonics  of  v and  i^,  . for  various  values 

o firms; 

of  L . 
a 

7.  TESTR2:  Calling  program  for  fault  current  simulation. 

8.  UNCONT:  Finds  the  solution  for  the  uncontrolled  rectifier  bridge  case. 
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CONT 


C MAIN  PLOTTER  PROGRAM 

DIMENSION  FD<5)  »F(5f  1 ) r FF(5»5)  »MUS(50)  rRMSIFOO) » 

1XRMS(50) f ALIF<21 ) fTBETA(21 ) »TMU(21 > »FF1 <4f<>  »F1 (4f 1 ) » 
1FDK4)  »ILS(50)  »BETAS<50)  »IFS(50)  »APIL(50)  »APBETA(50>  r 
1APMU(50) »APIF(50) »AHAR6(50) » AHAR12(50) »AHAR18(50) r 
1AHAR24(50) fAHAR30(50) »APRMS(50) fHAR6<50) »HAR12(50) » 
1HAR18(50) »HAR24(50) »HAR30(50) rAL0(21 ) »XHAR6(21 ) »XHAR12(50) » 
1XHAR18<21 > fXHAR24(21 ) r XHAR30 ( 21 ) r AIG( 60 > r AID(60) r AIF(60) f 
1AIK(60) .THETA(60) 

REAL*4  MOrMOO»MUfIF»ILrKlfLD»LF»MF»MFD»LAB»LO»MDrKQi- 
INNl fNP2r OMEGA f LA f MOP rMUFrKFOrKABrKODrlLSvlFSrMOSfKFfKD 
C INPUT  MACHINE  PARAMETERS 

WRITE<7»50) 

50  FORMAT( '0' f 2X» ' THIS  IS  THE  DATA  FOR  THE  PHASE  CONTROLLED 

1 BRIDGE  RECTIFIER' ) 

LF=0.12E  01 

LD=0.82E-07 

MF=0.79E-02 

MD=0.38E-05 

MFD=0.19E-03 

LAB=0.15E-03 

L0=0.3E-03 

C LO=LA  AND  LAB=MA 

,,  Kl  = 1.0 

} 0L=6760.0 

IL=1420.0 
0MEGA=2513.27 
FREQ=400.0 

KF= ( MF*LD-MD»MFD ) / ( LF*LD- < MFD ) **2 ) 
KD=(MD*LF-MF*MFD)/(LF*LD-(MFD)**2) 

1 KO=MD/LD 

I M00=MD*»2/LD 

MO=MOO 

DELTAF=MO+MOO 

I DELTAO=( 1 ,333*(L0+LAB) )-DELTAF 

' DIL=1420, 0/15.0 

KFO=MF/SQRT(LF»LO) 

KAB=LAB/LO 

KOD=MD/SQRT(LO*LD) 

I DL0=0.0 

IL=1420.0 
^ KKK=0 

f MF=KFO*SQRT(LF»LO) 

LAB=KAB*LO 

MD=KOD*SQRT<LO*LD) 

- M00=MD»«2/LD 

M0=M00 

DELTAF*MO+MOO 

DELTAO=( 1 .3333*<L0+LAB) )-DELTAF 
KF=(MF*LD-MD*MFD)/(LF*LD-<MFD)**2) 

CALL  NE WTON ( MU » BETA » I F f U > V » FREQ  $ DELTAO  r IL » VL » MOO  r K 1 1 
1 MF  r MO » OMEGA . DELTAF » ZETA ) 

IL=0.0 
IF=1  . li*IF 
DO  100  LLL»1»17 
KKK=0 
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CALL  PHACON<  IF  f HF  f PE TA  » hO » W . V » I L f MU » DELTAO  f OMEGA f VI.  f HOO ) 
APIL(LLL)  = Il. 

APPE1 A<LLL )=PETA*180.0/3. 1416 
APMU(LLL )=MU*1 80.0/3, 1416 

CALL  F S ( MU  f PE  I A f OMEGA  f 1 L » I F » MO  f DELT AO  f MF  » W » V f LO » LAB » 

1 DEL  T AF  » LLL  f AHAP6  » AHAR 1 2 » AHAR 1 8 > AHAR24  r AHAR30 ) 

CALL  RMS(BETAfLLL»MU»APRMSf IL»DELTAO»MF»MOfWfVFMOOr IFf 
IKF.DELTAF  ) 

IL-ILFDIL 
100  CONTINUE 

UIRI  rE(7f300)L0FlF 

300  FORMA r< '0' fSXf ' L0= ' » E 10 . 3 » 5X » ' IF* ' r E 1 0 , 3 ) 

UR I TEC  7f203) 

203  FORMAT < '0' f5Xr ' I L » IFRMS r 6TH f 1 2TH » 1 8TH f BETA f MU ' ) 

DO  201  KK=1f17 

WRITE(7f202)APIL(KK) fAPRMS(KK) fAHAR6(KK) fAHAR12(KK  ) f 
1AHAR18CKK) fAPBETA(KK) fAPMU(KK) 

201  CONTINUE 

202  FORMATC'  ' f 2X f F7 . 1 f 2X f F7 . 3 f 2X f E 10 . 3 f 2X f E 1 0 . 3 f 2X f 
1E10.3f2XfF7. 1 f2XFF5.2) 

STOP 

END 
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MACH? 


HI MENS ION  FD(5) »F<5f 1 ) » FF ( 5 f 5 ) f MUS ( 50 ) »RMSIF(50) » 

1XRMS(50) rALIF(21 ) »TBETA(21 ) »TMU(21 )»FFl(4f4)»Fl(4rl)» 

IFDl  <4)  f ILS(50)  riiETAS(50) » IFS(50)  »APIL(50)  »APBETA(50) » 
1APMU(50) »APIF(50) f AHAR6(50) »AHAR12(50) rAHAR18(50) » 
1AHAR24(50) fAHAR30(50) »APRMS(50) rHAR6(50) »HAR12<50) » 
1HAR18(50) rHAR24(50) »HAR30(50) fAL0(21 ) »XHAR6<21 ) »XHAR12<50) » 
1XHAR18(21 ) »XHAR24<21 ) »XHAR30(21 ) f AIQ(AO) »AID(&0) »AIF(60) f 
1AIK(60) f THETA(60) 

REAl.*4  MOfMOOfMUf  IFf  ILrKl  »LD»LF»MFrMFD»LAB»LOfMDfKQf 
INNl fNP2»0MEGA»LA»MUP»MUF»KF0fKAB»K0D» ILS» IFSfMUSt KF f KD 
C INPUT  MACHINE  PARAMETERS 

LF=0.12E  01 
LD=0.82E~07 
MF=0.79E-02 
MEi=0.38E-05 
MFD=0. 19E-03 
LAB=0. 15E-03 
L0=0.3E-03 

C LO=LA  AND  LAB=MA 

Kl=l ,0 
UL=6760.0 
IL=1420.0 
0MEGA=2513.27 
FREQ=400.0 

KF=(MF*LD-MD*MFD)/(LF*LD-<MFD)**2) 

KD=(MD*LF-MF*MFD)/(LF*LD-(MFD)*»2) 

KQ=MD/LD 

M00=MD*»2/LD 

MO=MOO 

DEl.TAF=MO+MOO 

DELTAO= ( 1 . 333* ( LO+LAB ) ) -DELTAF 
KFO=MF/SORT ( LF*LO ) 

KAB=LAB/LO 
KOri=MD/SQRT  ( LO*LD ) 

DL0=0.0 
H6M=9999.0 
DO  200  L=lf21 
LO=LO+DLO 
KKK=0 

MF=KFO*SQRT(LF*LO) 

LAB=KAB*LO 

MD=KOD*SQRT<LO*LD) 

M00*MD**2/l.D 

MO=MOO 

DELTAF=MO+MOO 

DELTAO=( 1 .3333*(L0+LAB) )-DELTAF 
KF= ( MF*LD-MD*MFD ) / ( LF*LD- < MFD ) **2 ) 

CALL  NEUTON<MUf BETA»IFfW»V»FREQfDELTAO»IL»VLfMOO»Kl » 

IMF  » MO » OMEGA » DELTAF » ZETA ) 

DL0=0.03E-03 
IF=1 . 1*IF 

CALL  PHACON( IF » MF r BETA f MO » W » IL » MU » DELTAO » OMEGA » OL f MOO ) 
CALL  FS(MUf  BETAfOMEGAf  ILf  IFtMOtDELTAOfMFfU»Vtl.OtLAB» 

1 DE L T AF  f L f AHAR6  f AHAR 1 2 f AHAR 1 8 f AHAR24  r AHAR30 ) 

CALL  RMS ( BETA f L f MU f APRMS f IL  f DELTAO f MF f MO f W f V f MOO f IF  f 
IKFf DELTAF) 

WR I TE ( 7 F 302 ) LO  f AHAR6  < L ) f APRMS ( L ) 
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302  FOKMAK '0' f2X. 'LA  = ' » E 1 0 . 3 f 5X » ' A TH  HARMONIC  OF  00  = 
lE10.3f5X. ' IF  RMS  = 'fElO.3) 

IF  ( H6M.lt.  AHARAd. ) ) GO  TO  200 
H6M=AHARA(L) 

BLA=L0 

BAPRMS=APRMS(L) 

200  CONTINUE 

WRI IE(7f 301 ) 

301  FORMAT ( '0' .5X» 'THIS  IS  THE  OPTIMUM  ARMATURE  INDUCTANCE') 

UR I TE ( 7 » 300 ) BLA » HAM » BAPRMS 

300  FQRMATC '0' r2Xr 'BEST  LA  = ' » E 10 . 3 f 5X f ' PEAK  6TH  HARM.  OF  00  ='f 
1E10.3f IXf 'PEAK  IF  RMS='fE10.3) 

LO=BLA 

MF=KFO*SQRT(LF»LO) 

LAB=KAB»LO 

MD=KOD*SQRT(LO*LD) 

l.AB=KAB*LO 

MD=KOD*SQRT(LO*LD) 

MOO=MH*MD/LD 

MO=MOO 

DELrAF=MO+MOO 

DELTAO=( 1 .3333*(L0+LAB) )-DELTAF 
KF= ( MF*LD-MD*MFD ) / ( LF*LD-MFD*MFD ) 

RETURN 

END 
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C MASTER  OPTOMIZATION  OF  L-C  FILTER  DESIGN 

REAL«4  K1 fK2f ILFHAXf IFfHUfLA»LADfLFfLO>MDfMFfMFDr 
1IL»MO»LWT»LEN»N1»LO 

21  CONTINUE 

WRITE(7. 1 ) 

WRITE(7f 104) 

104  FORMAT ( '0' »5X» 'ALL  INPUTS  HAVE  FORMAT  = F7.2  UNLESS 
1 OTHERWISE  SPECIFIED') 

1 FORMAT < '0' » IX » 'WRITE  THE  FOLLOWING  PARAMETERS  FOR  THE  FILTER') 
WRITE(7,2) 

2 FORMAT < '0' . IX » 'CAP.  ENERGY  DENSITY  < JOULES/LB.)  *') 

100  F0RMAT(F7.2) 

103  FORMAT( '0' f5X» 'F0RMAT=I2' ) 

101  F0RMAT(I2) 

READ(5»100)DC 

WRITE(7»4) 

4 FORMAT < '0' fix f 'CURRENT  DENSITY  FOR  L HIRE  (CIR  MIL/AMP)  =') 

READ<5f 100)CMA 
WRITE(7f 6) 

6 FORMAT( '0' f IXf 'MAX.  RMS  VALUE  OF  6TH  HARM.  OF  VO  (VOLTS)  »') 
READ<5f 100)V2 

WRITE(7f 7) 

7 FORMAT ( '0' f IX r 'ALLOWABLE  PEAK  FAULT  CURRENT  (AMPS)  =' ) 

READ(5f 100)ILFMAX 

IL=1420.0 
VI  =6760.0 
LF=1 .2 
LD=0.82E-07 
MFD=0. 19E-03 
0MEGA=2513.27 
no  50  IZ=lf2 

C IZ=1  IS  FOR  LA  NORMAL 

C IZ=2  IS  OR  LA  OPTIMUM 

IF(IZ.EQ.2)  GO  TO  51 
IF=275.0 
BETA=2. 12 
MU=0.307 
LA=0.300E-03 
LAB=0. 15E-03 

C LAB=MA 

MD=0.38E-05 
MF-0.79E-02 
W=144.0 
V=-138.0 
DELTA0=0.248E-03 
M0»0. 176E-03 
V1*0.120E  04 
WRITE(7f53)LA 

53  FORMAT ( '0' f2Xf 'THE  FOLLOWING  VALUES  ARE  BASED  ON  NORMAL 

ILA  =' fE10.3r2Xf 'H. ' ) 

GO  TO  52 

51  CONTINUE 

IL=1420.0 
IF=239.0 
nETA=2.28 
MU-0.604 
LA=0.720E-03 
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LAB=0.360E-03 
C LAP*MA 

MD=0.5B9E-05 
ME=0. 122E-01 
W=164.0 
y=-340.0 

DELTA0=0.595E-03 

M0=0.423E-03 

Ul=620.0 

WRITE(7»54)LA 

54  FORMAT ( '0' »2Xf 'THE  FOLLOWING  VALUES  ARE  BASED  ON 

lOPTIMUM  LA  =' »E10,3»2X» 'H. ' > 

52  URITE(7»300)IF»BETA»MU 

300  FORMAT ( ' ' r2X» ' IF=' »E10.3f5X» 'BETA=' fE10.3f5X» 'MU=' fElO.3) 

WRITE<7»305)IL»VL»V1 

305  FORMAT('  ' . 2X » ' I L= ' » E 1 0 . 3 r 5X » ' VL=> ' » E 10 . 3 r 5X » ' VI  = ' » E 10 . 3 ) 
K1=V2/V1 

C FIND  CONDUCTOR  AREA  IN  CIR  MILS  & SQ.  CM. 

CM=CMA»IL 

A16=(l+Ki )/(Kim36*0MEGA*»2) 

AM=CM*5.07E-10 

FW=AM 

F=0.7 

DW=5937.8 

A2=(VL*FW)**2*(K1+1 )/<DC*51 . 1E-07*K1«36*(0MEGA»*2)*F*F ) 
A3=3*3. 141A*F*DW 
A=(5*A2/(3*A3) )**0. 125 
N=A*»2»F/FW 

L0=(25.5E-07)»( (F/FW)**2)»(A)**5 
C0=A16/L0 

RES=N*3*3. 1416*A*(2.83E-08)/FW 
URITE(7r55)L0.C;0»RES 

55  FORMAT( '0' » IX. 'OPT. VALUES  BEFORE  FAULT  TEST 

1 ARE  L0=' .E10.3. IX. 'H. f CO* ' . ElO . 3 . 1 X . ' FD . ' . 1 X . ' RES= ' .E 10 . 3 ) 
K = 1 

15  RES=N*3*3. 1416»A»(2.83E-08)/FW 

CALL  FAULT2( IF.BETA.MU.LA.LAB.LF.LD.MD.MF.MFD.W.V.IL. 

1 OMEGA . DELTAO . MO . VL . LO . ILFMAX . RES ) 

WRITE(7.345)IL 

345  FORMAT('  '.2X.'MAX  LOAD  CURRENT  FROM  FAULT  ='.E10.3) 

IF( IL.LT. ILFMAX)  GO  TO  14 

K=K+1 

L0=L0*1 . 1 

C0=A16/L0 

A=( <L0*(FW/F)»*2)/25.5E-07)**0.2 
N=(A**2)*r/FW 
IF(K.GT.IOO)  GO  TO  14 
IL*1420.0 

IF(K.GT.2)  GO  TO  15 
UR1TE(7. 16) 

16  FORMAT( '0' . IX. 'FAULT  CURRENT  TOO  LARGE.  LO  INCREASED') 

GO  TO  15 
14  CONTINUE 

LUT-A3»A«*3 

CUT-A2/(A*«5) 

UT*LWT+CWT 

WRITE(7.17)L0.C0.IL.RES 
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FORMAT < 'O' »1X» 'LO=' »E10.3f5X» ' CO» ' » ElO . 3 »5X r ' ILF= ' »E10.3rlXr ' 

1 RES='fE10,3> 

WRITE<7»18>LMTfCWTfW7 

FORMAT < 'O'  fix » 'LWT=' »E10.3»5X, 'CWT=' ,E10.3f 5Xr 'TOTAL  WT  =' rElO.3) 
RAD=2«A 

URITE(7f 19>NfRA0f A 

FORMAT< 'O' f IXf 'NO.  TURNS  »'fI4»5Xf'L  RADIUS  ='r 
lE10.3r 'M' »5Xf 'L  LENGTH  - ' »E10. 3» 'M' ) 

CONTINUE 
gRITE(7f 20) 

FORMAT( '0' flXf 'WRITE  *0*  TO  ENDf  OR  M*  FOR  ANOTHER  RUN') 

WRITE(7f 103) 

REAIKSf  101  )KEY 
IF<KEY.GT.O)  GO  TO  21 
STOP 
END 
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C MAIN  PLOTTER  PROGRAM 

DIMENSION  FD(S) »F(5»1 ) fFF(5»5) »MUS(50) »RMSIF(50) f 
1XRMS(50>  rALIF<21 ) f TBETA(21 ) rTMU(21 ) ^ FFI < 4 f 4 ) »F1 ( 4 » 1 ) f 
IFDl <4) » ILS(50) » BETAS (50) f IFS(50) r APIL(50) »APBETA(50) f 
1APMU(50) »APIF(50) .AHAR6(50) .AHAR12(50) »AHAR18(50) » 
1AHAR24(50) »AHAR30(50) r APRMS ( 50 ) f HAR6 ( 50 ) r HARl 2 ( 50 ) » 
1HAR18(50) »HAR24(50) fHAR30(50) f AL0(21 ) »XHAR6(21 ) »XHAR12(50) r 
1XHAR18<21 ) fXHAR24(21 ) »XHAR30(21 ) rAIQ(60) f AI D < 60 ) » AIF ( 60 ) » 
1AIK(60) »THETA(60) 

RE AL»4  MO » MOO . MU » I F » I L » M f LD  * LF  f MF  f MFD » LAB  f LO » MD  f KCI . 

INNl »NP2»0MEGAfLA»MUP.MUFfKF0»KAB»K0D» ILSf IFSfMUSfKFfKD 
C INPUT  MACHINE  PARAMETERS 

LF=0.I2E  01 
LD=0.82E-07 
MF=0.79E-02 
MD=0.38E-05 
MFD=0. 19E-03 
LAB=0. 15E-03 
L0=0.3E-03 

C LO=LA  AND  LAB=MA 

Kl=l .0 
0L=6760.0 
IL=1420.0 
0MEGA=2513.27 
FREQ=400.0 

KF=(MF»LD-MD»MFD>/(LF»LD-(MFD)»»2) 

KD=(MD*LF-MF»MFD)/(LF*LD-(MFD)**2) 

KQ=MD/LD 

M00»MD**2/LD 

M0»M00 

DELTAF*M0+M00 

PELTA0=( 1 .333»(L0FLAB) )-DELTAF 
DIL=2130. 0/50,0 

CALL  NEWT0N(MU»BETA»IF»U.V»FRE0»DELTA0»IL»UL»M00»K1» 

1 MF  » MO  f OMEGA » DELTAF » ZET A ) 

CALL  LITLI(DELTAO»MF»IFfBETA»MO»Wf 
IMOOf V r IL r DELTAF »MU» AID» AIQ » AIFf AIK » THETA »KQ»KFfKD) 
KFO*MF/SQRT(LF*LO) 

KAB=LAB/LO 

KOD=MD/SQRT(LO*LD) 

UL0=*0.0 
1L»1420.0 
DO  200  L=»l  f 21 
L0»L0+DL0 
KKK=0 

MF»KF0*SQRT(LF*L0) 

LAB»KAB«L0 
MD<=KOO»SQRT  ( L0«LD ) 

M00»MD**2/LD 

M0=M00 

DELTAF=M0+M00 

DELTAO»( 1 ,S333*(L0+LAB) )-DELTAF 
KF=»  ( MF*LD-MD*MFD ) / < LF*LD-  < MFD  ) **2  > 

CALL  NE«TON<MU»BETA» IF»WfV»FREO»DELTAOf ILfOL»MOO»Kl r 
1 MF  f MO » OMEGA  r DELTAF » ZET A ) 

IF»1.1*1F 

CALL  PHACON(IF»MF»BETA»MO»W»V,ILfMUfDELTAO»OMEGA.ULrMOO> 
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ALO(L)=LO 

CALL  FS(MUfBETArOMEGAf ILrIFfHOfDELTAOrMFrUrVfLOrLAbf 
1 DEL  T AF  f L f XHAR6  r XHAR 1 2 r XHAR 1 8 1 XHAR24 » XHAR30 ) 

CALL  RHS(EiETArLfMUrXRHSf  IL»DELTAOfMF»MOrU»GfHOOf  IFt 
IKFfDELTAF) 

ALIF(L)=»IF 

TBETA(L)*BETA*180./3. 1416 
TMU(L)*MU»180./3.1416 
DL0=0.03E-03 
200  CONTINUE 

WRITE(7»201 ) 

DO  202  L=lf21 

URI TE(7>203>AL0(L>  »XRHS(L) fXHAR6(L) fXHAR12(L) »XHAR18(L) f 
lTBETA(L)fTMU(L) 

02  CONTINUE 

01  FORMAT<'  '»5X»'L0'»10X»'IF'rl0X»'6TH'»10Xf'12TH'fl0X»'18TH'f 

112X» 'BETA' f 12Xf 'MU' ) 

03  FORMAT('  '»lXf7E13.4) 

L0=0.3E-03 

LAB=0,15E-03 

MD=0.38E-05 

MOO*MD*MD/LD 

MO=MOO 

MF*0.79E-02 

KF= ( MF*LD-MD*MFD ) / ( LF*LD-MFD*MFD ) 

DELTAF=M0+M00 

tiELTAO*  ( 1 . 3333*  ( LO+LAB ) ) -DELTAF 
IL=0,0 

URITE(7f2800) 

DO  100  LLL*1»50 
KKK=0 

CALL  NEUTON<MU>BETA»IFfU»UfFREQ»DELTAO>ILfULfMOO»Klf 
IMF  f MO » OMEGA » DELTAF » ZETA ) 

ILS(LLL)»IL 

BETASCLLD-BETAtlOO. 0/3. 1416 
MUS(LLL)=MU*180. 0/3, 1416 
IFS<LLL)=IF 

800  FORMAT('  ' » 4X f ' IL ' » 8X f ' BETA' » 7X » ' MU' » 8X » ' IF ' ) 

801  FORMAT ( ' 0 ' » 2X f F7 . 1 »4X rF6 . 2»4X »F6. 3» 4X »F6 , 1 » 4X f F6 , 1 f 4X »F6. 1 ) 
CALL  FS(MUfBETAfOMEOAf ILFlFfMOfDELTAOrMFFUFUFLOrLABF 

1 DEL  T AF  F LLL  F HAR6  f HAR 1 2 f HAR 1 8 f HAR24  f HAR30  > 

CALL  RMS(BETAfLLLfMUfRMSIFfILfDELTAOfMFfMOfWfVfMOOfIFf 
IKFf DELTAF) 

IF=1 . 1*IF 

CALL  PHACON( IF f MF f BETA f MOfW f Of IL fMU f DELTAOf OMEGA fVL fMOO) 

APIL(LLL)=IL 

APBETA<LLL)=BETA 

APMU(LLL)=MU 

APIF(LLL)=IF 

CALL  FS(MUfBETAfOMEGAf ILfIFfMOfDELTAOfMFfWfUfLOfLABf 
IDELTAFfLLLf  AMAR6  f AHAR 1 2 f AHAR 1 8 r AH AR24  f AH AR30 ) 

CALL  RMS ( BET A F LLL  f MU f APRMS f I L f DELTAO f HF  f MO f W f U f MOO f IF f 
IKFfDELTAF) 

IL=IL+DIL 

IF(LLL.EQ.l)  IL»35.0 
100  CONTINUE 

DO  101  LLL=1f50 
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'2802 

102 

2803 

103 

104 

105 
1 10 

107 

106 


WRITE (7.2801 ) 1 LS ( LLL) . BETAS ( LLL ) .MUSCLLL ) . IFS(LLL ) 

CONTINUE 
DO  102  LLL-1.50 

WRITE(7.2802)LLL.HAR6(LLL> .HAR12<LLL) .HAR18(LLL) .HAR24(LLL) » 
1HAR30(LLL) 

FORMAT('  ' rlX»I3.5E13.4) 

CONTINUE 

DO  103  LLL=1.50 

WR I TE ( 7 . 2803 ) LLL . RMS I F ( LLL ) 

FORMAK'  ' . IX.  I4.L20.5) 

CONTINUE 
DO  104  L=1.50 

WR I TE( 7.2802) L. AHAR6(L > .AHAR12(L) .AHAR18(L  > . AHAR24 ( L ) . AHAR30 ( L ) 

CONTINUE 

DO  105  L=1.50 

WRITE (7. 2003 )L.APRMS(L) 

CONTINUE 
WRITE(7. 110) 

FORMAT( '0' .IX. 'IL.BETA.MU.IF— WITH  PHASE  CONT ' ) 

DO  106  L=1.50 

WRITE(7.107)APIL(L) . APBETA ( L ) . APMU ( L ) . APIF ( L ) 

FORMAT( '0' .1X.4E13.4) 

CONTINUE 

STOP 

END 
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C TEST  PROGRAM  FOR  NEWTON  SENSITIVITY  ANALYSIS 

DIMENSION  FD(5) rF(5»l ) fFFCSfS) f AIF ( 1 1 ) » ABETA( 1 1 ) rAMU< 11) r 
1AU< 1 1 ) r AV( 1 1 > f Al( 1 1 ) » AMO( 1 1 > > AMOO( 1 1 ) » ALAMD< 11 > 

REAl.*4  MOrMOOfMU»IFf  IL»KlfLD»LF»HFrMFD»LABfLO»MDfKQrlLFMAXf 
ILAfKFOfKABrKODfKFfKDrlKrMArlFTf IQfLAMQrLAMDf lLO»ILC»IKOf IKC 
DO  34  Il=l»7 
LF=0.12E  01 
LD=0.82E-07 
MF=0.79E-02 
MD=0.38E-05 
MFD=0.19E-03 
LAB=0.15E-03 
1.0=0. 3E-03 

C LO=LA  AND  LAB=MA 

CAF=MF/SORT ( LO*LF ) 

CAB=LAB/LO 

CAD=MD/SQRT(LO»LD) 

CFD=MFD/SORT  < LF*LD ) 

DO  35  12=1 » 11 
IF(I1.GT.6)  GO  TO  36 
IFdl.GT.S)  GO  TO  37 
IFdl.GT.A)  60  TO  38 
IFdl.GT.3)  GO  TO  39 
IFdl.GT.2)  GO  TO  40 
IFdl.GT.l)  GO  TO  41 
LF=0.6+0.12*d2-l> 

MF=CAF*SQRT  < LO*LF ) 

MFD=CFD*SQRT ( LF*LD  > 

A1(I2)=LF 
GO  TO  42 

41  LD=0.41E-07+<0.082E-07)*(I2-1) 

MD=CAD«SQRT ( LO«LO ) 

MFD=CFD*SQRT ( LF*LD ) 

AKI2)=LD 

GO  TO  42 

40  MF=0.395E-02+(0.079E-02)*(I2-l ) 

A1 ( I2)=MF 
GO  TO  42 

39  MFD=0 . 095E-03+ ( 0 . 01 9E-03 ) * ( 12-1 ) 

Ald2)»MFD 
GO  TO  42 

38  MD=0.19E-05+(0.038E-05)*(I2-l) 

A1 ( I2)*MD 
GO  TO  42 

37  LAB=0.075E-03+(0.015E-03)*(I2-l) 

A1 ( I2)=LAB 
GO  ro  42 

36  L0=0. 15E-03T<0.03E-03)*< 12-1 ) 

MF =CAF*S0RT ( LO*LF ) 

md=i:ad«sqrt(LO*ld> 

I AB=CAB«LO 
A1  ( I2)*l.O 

42  CONTINUE 
M«1  .0 
VL-6760.0 
IL-1420.0 
nMEOA-2513.27 
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FREQ=*400.0 

KF  = (HF*LD-MD*MrLi)/(LF*LD-MFD»MFD) 
KD=(MD*LF-MF*MFD)/<LF»LD-MFD*MFD) 

KQ=MD/LD 

M00=MD*»2/LD 

AM00<I2)=M00 

AMO ( 12 ) =KF*MF+KD*MD 

ALAMD( I2>=AM0( I2)-M00 

MO=MOO 

riELTAF=MO  + MOO 

DELTAO=( 1 .3333*<LO+LAB) )-DELTAF 
I CONTINUE 

CALL  NEW30N(MUfBETA»IF»W»0»FREQ»DELTA0r  ILfVl.fMOO.Kl  » 
IMF » MO f OMEGA »DELTAF.ZETA»K) 

IF(K.LT.200)  GO  TO  80 
J AIF(I2)=0.0 

ABETA( I2)=0.0 
AMU( I2)=0.0 
AW( I2)=0.0 
A0( I2)=0.0 
GO  TO  35 
) CONTINUE 

AIF(I2)=IF 
ABETA< I2)=BETA 
AMU( I2)=MU 
AU( I2)=W 
A0( IZ)’^ 

5 CONTINUE 

IF(I1.GT.6)  GO  TO  46 
IFdl.GT.S)  GO  TO  47 
IF<I1.0T.4)  GO  TO  48 
IF(I1.GT.3)  GO  TO  49 
IF( II .GT.2)  GO  TO  50 
IF(Il.GT.l)  GO  TO  51 
WRITE(7» 100) 

)0  FORMAT( '0' »5X» 'LF  VARIATION' ) 

GO  TO  52 

L WRITE(7»101) 

>1  FORMAT< '0' »5Xf 'LD  VARIATION' ) 

GO  TO  52 

) WRITE(7rl02) 

)2  FORMAT< '0' »5X» 'MF  VARIATION') 

GO  TO  52 

f URITE(7fl03) 

)3  FORMAT< '0' fSXf 'MFD  VARIATION' ) 

GO  TO  52 

J WRITE(7rl04) 

)4  FORHAT( '0' »5X» 'MD  VARIA1 ION' ) 

GO  TO  52 

' URITE(7»105) 

)5  FORMAT( '0' »5Xf 'LAB  VARIATION' ) 

GO  TO  52 

i WRITE<7»106) 

)6  FORMAT< '0' »5X» 'LO  VARIATION' ) 

I CONTINUE 

URITE(7fl20) 

20  FORMAT ( '0' »7X» 'PARAMETER ' 0X »' IF '. 1 IX »' BETA '» 1 IX t 'MU' » 


125 


110 

55 

34 


113X»'U'rl3Xf'V'»llX»' MOO ' f 1 2X » ' MO ' » 1 1 X » ' LAMD ' ) 

DO  55  L»l»ll 

URITE(7> 110)A1 (L) f AIF(L) f ADETA<L>  f AMU(L) >AU(L) f AO(L>  f 
lAMOO(L) rAMO(L) rALAMD(L) 

FORMAT ('  'f2Xr9E14.4) 

CONTINUE 

CONTINUE 

STOP 

END 
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C PROGRAM  TABLE. FOR 

1:  MAIN  PLOTTER  PROGRAM 

DIMENSION  FD<5) rP<5»I ) » FF ( 5 f 5 ) r MUS < 50 ) » RMSIF ( 50 ) » 

1XRMS(50) f ALIF(21 ) »TBETA(21 ) »TMU(21 )rFFl(4»4)»Fl(4»l)» 

• IFDl (4) f ILS(50) »BETAS<50) f IFS(50) rAPIL(50) r APBETA(50) r 

1APMU<50) »APIF(50) r AHAR6(50) »AHAR12(50) »AHAR18(50) r 
1AHAR24(50) rAHAR30(50) r APRMS(50) fHAR6(50) ,HAR12(50) r 
1HAR18(50) »HAR24<50) .HAR30(50) f AL0(21 ) »XHAR6(21 ) rXHAR12<50) r 
1XHAR18(21 ) »XHAR24(21 ) .XHAR30(21 ) » AI 0 ( 60 ) , AI D ( 60 ) r AIF ( 60 ) . 
1AIK(60) »THETA(60) 

REAL*4  MO»MOO»MUf IF» ILrM  »LD»LFrMFfMFD»LAB»LO»MDfKQ» 

INNl »NP2»0MEGAf LArMUP»MUF»KF  0»KABfK0D» ILS» IFSfMUSfKF  rKD 
C INPUT  MACHINE  PARAMETERS 

LF=0.12E  01 
LD=0.82E-07 
MF=0.79E-02 
MD=0.38E-05 
MFD=0. 19E-03 
LAB=0.15E-03 
L0=0.3E-03 
Kl  = l .0 
0L=6760.0 
IL=1420.0 
0MEGA=2513.27 
FRE0=400.0 

KF= ( MF»LD-MD»MFD ) / ( LF*LD- ( MFD ) **2 ) 
KD=(MD»LF-MF*MFD)/(LF»LD-(MFD)**2) 

KQ=MD/LD 

M00=MD»*2/LD 

MO=MOO 

DELTAF=MO+MOO 

DELTAO=( 1 .333*(L0+LAB) )-DELTAF 
DIL=1420. 0/15.0 
KFO^MF /SORT  (LF  *1.0) 

KAB=LAB/LO 

KOD=MD/SQRT(LO*LD) 

DL0=0.0 

DO  200  L*l»21 

IL»1420.0 

LO=LO+DLO 

KKK*0 

MF=KF0*SQRT(LF*L0) 

LAB»KAB*L0 

MD*KOD*SQRT(LO*LD) 

M00=MD**2/LD 

MO-MOO 

DELTAF-MO+MOO 

DELTAO=( 1 .3333*(L0+LAB) )-DELTAF 
KF-(MF*LD-MD*MFD)/(LF*LD-(MFD>**2) 

CALL  NEWTON(MU»BETA»IFfWfV»FREQ»DELTAO. ILfVLfMOOfKl f 
1 MF  » MO » OMEGh • DEI  T AF  t ZETA ) 

DL0-0.03E-03 

IL-O.O 

IF-1 . 1*IF 

DO  100  LLL*1»17 

KKK-0 

CALL  PHACON( IFfMF»BETA»MO»W»V.IL»MU»DELTAO»OMEGA»VL»MOO) 
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APIL(LLL)=IL 

APBETA(LLL)=BETA» 180. 0/3. 1416 
APHU(LLL)=MU*180.0/3. 1416 

CALL  FS(MU.BETA.OMEGA»ILf IF»M0»DELTA0»MF»W»V»L0fLAB» 

1 DEL T AF  f LLL  f AHAR6 » AHAR 1 2 » AHAR 1 8 » AHAR24  r AHAR30 ) 

CALL  RMS ( BE TA  r LLL » MU » APRMS r I L f DELTAO » HF » MO » W r U » MOO » I F f 
IKF.DELTAF) 

IL=IL+DIL 
100  CONTINUE 

WRITE(7r300)L0» IF 

300  FORMAT ( '0' »5X. 'L0=' »E10.3r5X» ' IF*' rElO.3) 

URITE(7,203) 

203  FORMAT( '0' f5X» ' IL . IFRMS.6TH» 12TH» IBTHrBETAf MU' ) 

DO  201  KK=1»17 

URITE(7.202)APIL(KK) . APRMS ( KK )» AHAR6 ( KK ) »AHAR12(KK) t 
1AHAR18(KK) »APBETA(KK) »APMU(KK) 

01  CONTINUE 
00  CONTINUE 

02  FORMAT!'  ' r 2X . F7 . 1 . 2X f F7 . 3 r 2X » El 0 . 3 . 2X f ElO . 3 » 2X . 
lEl 0. 3. 2X rF7. 1 f2X.F5. 2) 

STOP 

END 
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u") 


Tt;siK? 


C lEST  PROGRAM  FOR  FAULI 

DIMENSION  FD(5f  1 ) fF<5»  1 ) »FF(5»5)  f A(4»n)  ,P(5f5)  fRH(5f  1 ) »CIl.  (300)  , 
1CIK< 300) »CWT(300) 

REAL*4  MO f MOO, MU » IF , I L , M , L D , LF r MF r MFD , I AB , LO , MU , KO r ILFMAX , 

ILA,  KFO,KAD,KOn,KF,Krtf  IKfMAr  IFF,  IQfl  AMO,  1.  AMD,  II  0 r 1 l.C  , I KO » IKf. 
LF=0.12E  01 
LD=0 .B2E-07 
MF=0.79E-02 
MD=0.38E-05 
MFD=0. 19E-03 
LAB=0. 15E-03 
L0=0.3E-03 
l.A=LO 

ILFMAX=0.1E  06 
URITE(7r200) 

200  FORMAT ( '0' ,2Xf 'TYPE  THE  OALUE  OF  RP  DESIRED') 

READ(5, 201 )RP 

201  FORMAT (E20. 10) 

K 1 = 1 . 0 

0L=6760.0 

IL=1420.0 

0MEGA=2513.27 

FREQ=400.0 

KF=(MF*LD-MD»MFD)/(LF»LD-MFD*MFD) 

KD=  ( MD*LF-  MF*MFD ) / ( LFXrLD-Mr  D*MFD ) 

KQ=MD/LD 

M00=MD»*2/LD 

M0=M00 

DELTAF=M0+M00 

DELTAO=( 1 .3333*(L0FLAB) )-DFLTAF 

CALL  NEWTON(MUf BETA, IF,W,0,FRFQrDELTA0, IL,VL,M00,K1 , 

IMF, MO, OMEGA »FEl TAF,ZEIA) 

00  URirE<7,299) 

99  FORMAT( '0' , IX, 'ENTER  OALUE  OF  LO  DESIRED FORMAT  E10.3') 

READ(5,298)L0 
298  FORMAT (ElO. 3) 

10  FORMAK'  ',1X,5E13.4) 

IL=1420.0 

CALL  FAULT2( IF , BETA , MU , LA , L AB , LF , I D , MD , MF , MFD , W , V , IL, 

1 OMEGA, DELTAO, MO, OL,LO, ILFMAX,RP) 

WRI TE( 7,300) IL 

300  F ORMAT( '0' ,3X, ' ILMAX  FROM  FAULT  =',E10.3) 

WRITE( 7,600) 

600  FORMATC '0' ,3X, 'TYPE  'O*  TO  END,  OR  *1*  FOR  ANOTHER  RUN') 

RE AD (5, 700) KEY 
700  FORMATC 14) 

IFCKEY.GT.O)  GO  TO  500 

STOP 

END 


129. 


UN CON T 


C MAIN  PLOTTER  PROGRAM 

HI  MENS  ION  rri<5)  »r  (5»  1 ) »FF(5»5)  »MUS(50)  fRHSIF(50)  f 
1XRMS(50) »ALIF(21 ) »TBETA(21 ) »TMU(21 )»FFl(4»4)fFl(4»l)» 

IFDl (4) f ILS<50) »BETAS<50) f IFS(50) » APIL ( 50 ) » APBETA < 50 ) » 
lAPHU(SO) >APIF(50) f AHAR6(50) r AHAR12C50) r AHAR18(50) f 
1AHAR24(50>  f AHAR30<50) f APRMS ( 50 ) » HAR6 ( 50 ) »HAR12(50) f 
1HAR18(50) »HAR24(50) »HAR30(50) »AL0(21 ) rXHAR6(21 ) » XHARl 2 ( 50 ) f 
1XMAR18(21  > fXMAR24(21  ) fXHAR30(21  ) fAIQ(60)  f AID ( 60 ) f AIF  ( <S0 ) f 
1AIK(60>  »THETA(60) 

REAL*4  MOfMOOfMUf IF  f IL f K 1 f LD f LF f MF f MFD f LABf LOf MD f KQ f 
INNl fNP2f OMEGA F LA fMUPfMUFfKFOfKABfKODfILSfIFSfMUSfKFfKD 
C INPUT  MACHINE  PARAMETERS 

WRITE (7f50) 

50  FORMAT( 'O' f2Xf ' THIS  IS  THE  DATA  FOR  THE  UNCONTROLLED 
1 BRIDGE  RECTIFIER'  ) 

LF=0.12E  01 

LD=0.82E-07 

MF=0.79E-02 

MD=0.38E-05 

MFD=0. 19E-03 

LAP=0. 15E-03 

L0-0.3E-03 

C LO=LA  AND  LAB=MA 

K 1 = 1 . 0 
0L=6760.0 
IL=1420.0 
0MEGA=2513.27 
FREQ=400.0 

KF=(MF*LD-MD*MFD)/(LF*LD-(MFD)**2) 

Mi=(MD*LF-MF*MFD>/<LF»LD-(MFD)**2) 

KQ=MD/LD 

M00=MD**2/LD 

MO=MOO 

DELTAF=MO+MOO 

DELTAO=< 1 ,333*(L0+LAB) )-DELTAF 
DIL=1420. 0/15.0 
KFO=MF/SQRT  < LF*LO ) 

KAB=LAB/LO 

K0D=MD/SQRT(L0*LD) 

HL0=0.0 
IL=1420.0 
KKK  = 0 

Mf =KF0*S0RT(LF*L0) 

LAB=KAB»LO 
MD=KOD*SQRT (LO»LD) 

M00=MD**2/LD 

MO=MOO 

DEL  I AF=MO+MOO 

DELrAO=( 1 .3333*(L0+LAB) )-DELTAF 
KF= ( MF*LD-MD»MFD ) / ( LF*LD- ( MFD ) **2 ) 

IL=0.0 

DO  100  LLL=1f17 
KKK=0 

CALL  NEUTON(MUfBETAf IFfUfOfFREOfDELTAOfILfOLfMOOfKI F 
IMF  f MO  F OMEGA  f DELTAF  f ZETA ) 

Af IL<1 LL )*IL 
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AFIF(LLL)=IF 

AFPE7A(LLL)=EttTA*180.0/3. 1416 
APMU<LLL)=MU* 180. 0/3, 1416 

CALL  FS<MU»BETA.8ML8Af ILf IF  rMOf DELTA0»MF»W»0,L0»LAB» 

1 DEI  T AF  F LLL  f AHAR6  f AMAR 1 2 f AMARl 8 f AHAR24  f AHAR30 ) 

CALI.  RMS(BETAfLI  LfHUfAPRMSf  IL  f DELTAO  f MF  f MO  f W f 0 f MOO  f IFf 
IKFfDELTAF) 

IL=IL+DIL 
100  CONTINUE 

WRI TE(7f300)L0 

jOO  format < 'O' f5Xf 'L0=' fE10.3) 

URITE( 7f203) 

203  FORMAT  ( 'O'  f5Xf  ' 1 1.  f I FRMB  f 6TM  f 1 2TH  f 1 8TH  f BET  A f MU  f IF  ' ) 

DO  201  KK=1f17 

URITE(7f202)APIL(KK) fAPRMS(KK) fAHAR6(KK) fAHAR12(KK)  f 
1AHAR18(KK)  fAPE«ETA(KK)  f APMU  ( KK  ) f API  F ( KK  ) 

201  CONTINUE 

202  F0RMAT('  ' f 2X f F7 . 1 f 2X f F7 . 3 f 2X f E 1 0 . 3 f 2X f E 1 0 . 3 f 2X f 
1E10.3f2XfF7. 1 f2XfF5.2f3XfF6.2) 

STOP 

END 


131. 


APPENDIX  IV.  SUBROUTINES 


The  following  subroutines  for  the  main  programs  of  Appendix  III 
are  listed  in  alphabetical  order.  Subroutines  GELG  and  ARSIN  are  not 
included.  GELG  is  a program  for  solving  simultaneous  equations  that 
is  part  of  the  IBM  Scientific  Subroutine  Package.  ARSIN  is  a series 
for  the  arcsin  function.  It  should  be  noted  that  the  notation  in  the 
programs  occasionally  varies  from  that  in  the  text: 


1.  rS:  finds  the  harmonics  of  v . 

— o 

2.  FA’JLT.? : Cilculates  the  fault  current. 

3.  JACOB : Calculates  the  Jacobian  matrix  for  the  uncontrolled 

rectifier  bridge. 

4.  JAC0B4 : Calculates  the  Jacobian  matrix  for  the  controlled 

rectifier  bridge. 

5.  LITLI:  Calculates  i,,  i , i-  and  i,  vs.  B. 

d q f k 

6.  NEWTON : Newton-Raphson  algorithm  for  the  uncontrolled 

rectifier  bridge. 

7.  NEW30N : Same  as  NEWTON  except  variable  K is  included  in 

argument  list  to  test  for  convergence.  Used  only 
with  SENSI3. 

8.  PHACON : Newton-Raphson  algorithm  for  the  controlled 

rectifier  bridge. 
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9.  RHS : Calculcites  right  hand  side  vector  for  NEWTON. 

10.  RHS4B4 : Calculates  right  hand  side  vector  for  PHACON. 

11.  FM.S : Find  rms  value  of  i^.. 

12.  TERMA : Performs  I'epetitive  calculation  for  FS. 
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SUBROUTINE  FSCMUf BETArOHEOAr ILf IFfMO»DELTAO»MFrWfUfLO»LABr 
1 DEL  T AF  f LLL  r HAR6  f HARl 2 f HAR 1 8 r HAR24  r HAR30 ) 

DIMENSION  CN<5>  »HAR6<50>  fHAR12(50>  »HAR18(50>  »HAR24(50>  fHAR30(50> 

REAL»4  MUrIL»IFfMO»MFfLO»LAB 

A=-OMEGA*(l .732*IF*MF+2.865*M0*W)*1.91 

B=2 . 86S»0MEGA*M0$U«1 .91 

C=-2 . B65*0MEGA*IL»M0*l .91 

DD= ( OMEGA/DELTAO ) » ( 1 . 5»M0-L0-LAB  > 

D=DD»( 1 . 155*MF*IF+1 .91*M0*W)*1 .91 

E=+l .91*M0*V*DD*1 .91 

F=0,955*IL»DELTAF*DD*1 .91 

DO  10  K=lf5 

B1=BETA+MU-1 .047 

B2=B1+1 .047 

N=K*6 

CALL  TERMA(N»Ar2.094fB2f ANl fBNl > 

AN>AN1 

BN«=BN1 

CALL  TERMA(N>Ar2.094fBlfANlrBNl> 

AN-AN-ANl 

BN=BN-BN1 

CALL  TERMA(NfBr0.524fB2f ANlfBNl) 

AN=ANTAN1 

BN=BN+BN1 

CALL  TERHA(NfBf0.524rBl» ANlfBNl) 

AN=AN-AN1 

BN=BN-BN1 

ANG=2. 094-BETA-MU 

CALL  TERMA(NfCtAN0fB2fANl»BNl) 

AN=AN+AN1 

BN=BN+BN1 

CALL  TERMA(NrCfANG»BlfANl»BNl> 

AN=AN-AN1 

BN=BN-BN1 

B1=BETA 

B2=BETA+MU 

CALL  TERMA(NfDf0.0fB2fANlfBNl) 

AN«ANTAN1 

BN«BN-fBNl 

CALL  TERMA(NfDfO.OfBlfANl»BNl) 

AN-AN-ANl 

BN»BN-BN1 

CALL  TERMA(NfEfl. 571 »B2f ANlfBNl) 

AN=cANTANl 

BN«BN+BN1 

CALL  TERMA(NfEfl.571fBlfANlfBNl) 

AN=AN-AN1 

BN»BN-BN1 

ANG=-BETA-MU 

CALL  TERMA(NfFfANGfB2fANlfBNl ) 

AN=ANTAN1 

BN=BN+BN1 

CALL  TERMA(NfFfANGfBlfANlfBNl) 

AN»AN-AN1 

BN*BN-BN1 

CN<K)-S0RT<AN**2+BN**2) 
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CONTINUE 

HARA(LLL) 

MAR12(LLl. 

HAR18<LLL 

HAR24(LLL 

HAR30(LLL 

RETURN 

END 


CN(1  ) 

=CN( 

=CN( 

=CN( 

=CN( 


n ▼ in 


SUBROUTINE  FAULT2( IFf BETArMUrLAf LABf LFrLD»HDfHF»HFOf UrU» ILr 
1 OMEGA » BELT AO » MO  » VL  » LO » I LFMAX  » RP ) 

C ROTOR  FLUX  LINKAGES  ASSUMED  CONSTANT 

DIMENSION  A(Af4>  fB(5fS)  fCILOOO)  rCIKOOO)  rCWT(300)  r 
lARH<4f l)rBRH(Srl) 

REAL*4  LO»KDfKF»KQ» IK» IFfMUrLA»LAB»LF»LD»MDrMFfMFDf IL» IFTr 
lIDrIQrLAMQfLAMFrLAMDr ILOf ILCf IKOrlKCrHOrILFMAX 
KK=0 
JJ=0 

C SPECIFY  INITIAL  CONDITIONS  FOR  CONDUCTION  PERIOD 

WT=BETA+MU-1 .047 
C L0=0.1 

KD= ( MD*LF-MF*MFD ) / < LF*LD-MFD*HFD ) 
KF=(MF*LD-MD*MFD)/(LF*LD-MFD*MFD) 

KQ=MD/LD 

IFT=IF+1 .732»KF*( (3/3. 1416 >*(W+IL* 

1C0S(BETA+MU) ) - ( IL*COS ( WT+0 . 524 ) ) > 

ID=( IFT-IF)*KD/KF 

IQ=1 .732*KQ*( (3/3. 1416)*(V-IL*8IN(BETA+MU) ) 
l+IL*SIN(WT+0.524) ) 

F1*1/(LF*LD-MFD*MFD) 

LAMQ=- 1 . 732»MD*IL*SIN( UT+0 . 524  > f LD» IQ 
LAMF=1 .732»IL»MF*C0S( WT+0. 524 >+LF*IFT+MFD»ID 
LAMD=1 .732*IL*MD*C0S(WT+0.524)+MFD*IFT+LD*ID 
SF*F 1 ♦ ( LAMF*LD-LAMD*MFD ) 

SD=F 1 * ( LAMD*LF-LAMF*MFD ) 

SQ=»LAMQ/LD 

A(1»1)*L0+2*(LA+LAB) 

A(2»2)=l .0 
A(3f3)=»l  .0 
A(4f4)=l .0 
A(2»3)=0.0 
A(2f 4)=0.0 
A(3»2)=0.0 
A(3r4)=0.0 
A(4»2)=0.0 
A(4»3)=0.0 
B(1»1)=L0+2*(LA+LAB) 

B(2»2)=l .0 

B(3r3)=l .0 

B(4»4)=l .0 

B(2f3)=0.0 

B(2»4>=0.0 

B(3f2>=0.0 

B(3r4)>0.0 

B(4f2):>0.0 

B(4>3)30.0 

B(5f 1 )=-LA-LAB 

B(5f5)=2*(LA+LAB) 

B( 1 »5)=-LA-LAB 

C GO  INTO  CONDUCTION  DO  LOOP 

DWT=>=3. 1416/150 
DT=DWT/2513.3 
WT=WT+0.524 
)UMP*0 

110  WT=WT-1 .047* JUMP 
JUMP=JUMP+1 
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IF(  JUMF^.GK.4)  GO  TO  99 

RA=0.0164 

DO  1 K=1.3000 

KEY=0 

ILO=IL 

•12  KEY=KEY  + 1 

IF(KEY.GT.2)  GO  TO  71 
C FIND  DY/DX  AT  0 

A(  1»2)  = 1 .732*MF»COS(UT) 

A<  1 »3)  = 1 .732*MD*COS(WT) 

A(  1 f4)=-l ,732*MD*SIN(WT) 

A(2» 1 )=1 ,732*KF*COS<WT) 

A<3» 1 )=+l .732*KD»COS(UT ) 

A(4»  1 )=-l ,732»KO*SIN(WT) 

71  CONTINUE 

ARH( 1 f 1 )=-(RP+2»RA)*IL+l .732*MF»0MEGA*SIN(UT)* 
lIFT+1 .732»MD*OhEGA*SIN(WT)*ID 
141 .732»MD»0MEGA»C0S(UT)*IQ 
ARH<2rl )=1 .732*KF»0MEGAtSIN<UT)*IL 
ARH<3f 1 )=+l .732»KD*0MEGA»SIN(WT)*IL 
ARH(4» 1 )=1 ,732*K0*0MEGA*C0S(UT)*IL 

ARH( 1 » 1 )^( ARH( 1 » 1 )-A( 1 f 2)»ARH(2f 1 )-A< 1 f 3)*ARH(3f 1 ) 

1-A( 1 f 4 )*ARH<  4. 1 ))/(A(l»l)-A(l»2)*A(2»l) 

1-A< 1 »3)*A<3» 1 )-A( 1 f 4)*A(4» 1 ) ) 

ARM ( 2 » 1 > =ARH (2»1)-A(2»1) * ARM  < 1 » 1 ) 

ARH(3. 1 )=ARH(3» 1 )-A(3» 1 )*ARH( 1 »1 ) 

ARM ( 4 » 1 ) =ARH (4»1)-A<4rl) *ARH ( 1 r 1 ) 

£ ARH( 1 f 1 )=DIL/DT  AT  0 

IF(KEY.GT.l)  GO  TO  2 

UT=WT+DWT 
IL=IL4-ARH(  1 . 1 )*Df 
ILC=IL 

DILDT*ARH( 1 » 1 ) 

GO  TO  3 

2 CONTINUE 

IL=ILO+( (DILDT+ARH( 1.1) )*DT)/2 

3 CONTINUE 

IFT*SF-1 .732*IL»COS(WT)*KF 
ID=SD-1 . 732*IL*C0S(WT)*KD 
IQ*=5Q41 .732*KQ«IL*SIN(UT) 

IF<KEY.LT.2)  GO  TO  12 
IF(ABS(IL-ILC).LE.1.0)  GO  TO  10 
ILC=IL 

IF(KEY.0T.50)  GO  TO  997 
GO  TO  12 
10  CONTINUE 

KK»KK4-1 
CIL(KK)-IL 

IF(  IL.QT. ILFMAX)  GO  TO  997 
CIK<KK)-0.0 

IF(KK.0E.299)  GO  TO  997 
301  CONTINUE 

C TEST  FOR  END  OF  CONDUCTION  PERIOD 

VBCT-RA*IL+(LABfLA)*ARH(lrl )-1.732*HF*SIN<WT-0.524)*ARH(2.1) 

1-1 .732*0HE0A*C0S<WT-0.524)*(HF*IFT+MD*ID) 

1-1 ,732»MD*SIN(WT-0.524)*ARH(3rl )-l .732*MD*C0S(WT-0.524)*AFM(4»1 > 
1+1 .732*HD*0ME0A*SIN(WT-0.524)*IQ 
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C FAULT  IS  PRESENT  FOR  1 CONDUCTION  PERIOD  PLUS  1 

C COMMUTATION  PERIOD  + 1 COMMUTATION  PERIOD  WHERE  NEXT  SCRS 

C ARE  BLANKED— PROGRAM  ENDS  WHEN  READ  IL  IS  PAST 

IF((JUMP.GE.2>.AND.(IL.LT.CIL(KK-1))>  GO  TO  99 
IF< JUMP.GE.2)  GO  TO  1 
IF(VBCT.GT.O.O)  GO  TO  11 

I CONTINUE 
WR1TE(7»720) 

720  FORMAT('  ' r 5X f ' CONDUCTION  PERIOD  DOES  NOT  END') 

GO  TO  997 

II  CONTINUE 

C CALCULATE  COMMUTATION  INTERVAL 

IK=0.0 

DO  26  K=l»200 
KEY=0 
ILO=IL 
IKO=IK 

120  KEY=KEY+1 

IF(KEY.GT.2)  GO  TO  70 
C FIND  DY/DX  AT  0 

B< 1 f2)=l .732*MF»C0S(WT) 

B< 1 »3)=1 .732*MD*C0S(WT) 

B( 1 r4)=-l ,732»MD»SIN<WT) 

B(2»l )=1 ,732»KF»C0S(WT) 

B<2»b)=l .732*KF*SIN(WT-0.524> 

B(3»l )=+l .732*KD*C0S(WT) 

B(3f5)=+1 .732*KD*SIN(WT-0,524) 

B<^f 1)=-1.732*KQ*SIN(WT) 

B(4,5)=l .732*KQ*C0S(WT-0.524) 

B(5r2)=*l  .732*MF*SIN(WT-0.524) 

B(5f3)=l .732«MD«SIN(WT-0.S24) 

B(5»4)=l .732*MD*C0S(WT-0.524) 

70  CONTINUE 

BRH( 1 f 1 )=-(RPT2«RA)«ILTRA»IKf 1 .732»MF«0HEGA*SIN(WT) 

14IFTT1 .732«MD*0MEGA*SIN(WT)«ID 
1+1 .732*MD*0ME0A»C0S(WT)«IQ 
BRH(2fl )=1 .732*KF»0MEGA*(SIN<WT)*IL-IK*C0S(WT-0.524) ) 
BRH(3»1  )>=+l  .732«KD»0MEGA*(IL»SIN(WT)-IK«C0S(WT-0.S24) ) 

BRH ( 4 » 1 ) = 1 . 732*KQ*0MEGA* ( IL*COS ( WT ) +IK*SIN ( WT-0 . 524 ) ) 
BRH(5f 1 )°RA»IL-2«RA«IK-1 .732»MF«0ME0A«IFT* 

1C0S( WT-0 . 524 ) -1 . 732*MD»0MEGA*ID* 
lC0S(UT-0.524)+1.732»MD*0MEGA*ia* 
lSIN<UT-0.524) 

H=>B<lrl)-B(lf2)«B(2Fl>-B(l»3)»B(3fl)-B(l»4)»B(4»l) 

C=B( lf5)-B(lt2)*B(2»5)-B(lf3)*B(3f5)-B(l»4)*B<4f5) 

D=BRH( 1 »1 )-B( 1 »2)»BRH(2f 1 )-B( 1 » 3)*BRH( 3» 1 )-B< 1 »4)*BRH(4r 1 ) 
E»B(5f 1 )-B(5>2)«B(2r 1 )-B(5f3)«B(3rl )-B(5r4)«B(4rl ) 
F=B(5f5)-B(2»5)*B(5f2)-B(5f 3)*B(3»5)-B(5f4)*B(4»5) 
G=BRH(Srl)-B(5f2)«BRH(2f l)-B(5»3>«BRH(3f l)-B(5f4)»BRH(4rl> 
BRH( 1 »1 )=(D»F-C*G)/(H*F-C*E) 

BRH(5»1 )=<H»G-D*E)/(H*F-C*E) 

C BRH( 1 f 1 )=DIL/DT»  BRH(5r 1 )-DIK/DT  AT  0 
IF(KEY.GT.l)  GO  TO  20 
WT=UT+DWT 
IL*IL+BRH(1»1)*DT 
IK=IK+BRH(5f 1 )*DT 
ILC-IL 
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IKC=IK 

DILDT=BRH( 1 » 1 ) 

DIKDT=BRH(5» 1 ) 

GO  TO  30 
20  CONTINUE 

IL=ILO+( (DILDT+BRM( 1 » 1 ) )*DT)/2 
IK=IKO+( (DIKDT+BRH(5»1 ) )*DT)/2 
30  CONTINUE 

IFT=SF-1 .732*KF*( IL»C0S(UT)+IK*SIN(WT-0.524) ) 

ID=SH-1 .732*KD*( IL*C0S<WT)+IK*SIN(WT-0,524) ) 

IQ=S0+1 .732*Ka»(IL*SIN<WT)-IK»C0S(WT-0.524) ) 

IF(KEY.LT.2)  GO  TO  120 

IF( (ABS< IL-ILC) .LE. 1 .0) .AND. <ABS( IK-IKC) .LE. 1 .0) ) GO  TO  100 

ILC=IL 

IKC=IK 

IF(KEY.GT.SO)  GO  TO  997 
GO  TO  120 
100  CONTINUE 
KK=KK+1 
CIL(KK)»IL 
CIK(KK)»IK 
300  CONTINUE 

IF(IL.GT.ILFMAX)  00  TO  997 
C TEST  FOR  END  OF  COMMUTATION  PERIOD 

IF(IK.GE.IL)  GO  TO  110 
26  CONTINUE 

WRITE(7f721 ) 

■’21  FORMAT('  ' rSXf 'COMMUTATION  PERIOD  DOES  NOT  END') 

GJ  TO  997 
99  CONTINUE 

, DUM=»BETA+MU-1 .047 

997  DO  86  K»1»KK 

UUM=DUM+DMT 
CWT(K)=DUM 
86  CONTINUE 

DO  144  L*lrKK 

C WRITE ( 7» 145)CIL (L) »CIK(L) »CWT(L) 

145  FORMAT('  ' » 2X  f ' IL<= ' f E 10 . 3f  5X  f ' IK=  ' »E10 , 3 » 5X  f ' WT» ' r ElO . 3 ) 

144  CONTINUE 

RETURN 
END 
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SUBROUTINE  JACOB(HF » BETA» IF » Wf V» IL »MO»MU» DELTAOrOMEGAr FF » A»B»C) 
DIHENSION  FF(5f5) 

REAL«4  MOFHOOfMUrHFr ILfIFfKI f OMEGA 
PAIF=1 .155*MF*SIN(BETA) 

PBIF=1 . 155*MF 

PABETA=1 . 155*IF*MF»COS(BETA)+l .91*H0*<M*C0S(BETA)-V*SIN(BETA) ) 
PBBETA=-1 .91»M0*IL*SIN<BETA+MU) 

PCBETA=-1 .91*M0*IL*C0S(BETA+MU) 

PAMU=-1 .91*M0*IL*C0S(MU) 

PUMU=0.0 

PUMU=0.0 

PBMU=-1 ,91*M0*IL*SIN(BETA+MU> 

PCMU=-1 .91*M0*IL*C0S(BETA+MU) 

PAW=1 ,91*MO*SIN(BETA) 

PBU=1 .91)»M0 

PAV=1 .91»M0*C0S(BETA) 

PCU=1 .91*M0 

FF ( 1 » 1 ) = ( COS ( BETA+MU  > -COS ( BETA ) ) *PAIF 
1F0.25*(2*MU-SIN(2*(BETA+MU) )+SIN<2*BETA) )*PBIF 
FF(1 »2)=(C0S(BETA+MU)-C0S(BETA> )*PABETA 
1 -A* (SIN (BETA+MU) -SIN( BETA >) -0.25* (2*MU 
1-SIN(2*(BETA+MU) ) +SIN( 2*BETA ) )*PBBETA 
1+0.5*B*(-C0S(2*(BETA+MU) )+C0S(2*BETA) ) 

1+0. 5*( (SIN(BETA+MU) )**2-(SIN(BETA) )**2) 

1*PCBETA+C* (SIN( BETA+MU )*COS( BETA+MU )-SIN( BETA )*COS( BETA ) ) 

FF ( 1 » 3 ) =PAMU* ( COS ( BETA+MU ) -COS ( BETA  > ) -A*SIN ( BETA+MU ) 

1 +0. 25* (2*MU-SIN( 2* (BETA+MU) )+SIN(2*BETA) )* 
lPBMU+(B/2)*( 1-C0S(2*( BETA+MU) ) ) +0. 5* ( (SIN( BETA+MU) )**2 
1-(SIN(BETA) )**2)*PCMU+C*(SIN(BETA+MU)*C0S(BETA+MU) ) 

FF( 1 f4)=PAV*(C0S(BETA+MU)-C0S(BETA) )+0.5*( (SIN(BETA+MU) )**2 
1-(SIN(BETA) )**2)*PCV 

FF  ( 1 . 5 ) =DEl.TAO+PAW*  ( COS ( BETA+MU ) -COS ( BETA ) ) 
1+0.25*(2*MU-SIN(2*(BETA+MU) )+SIN(2*BETA) )*PBW 
FF(2»1 )=PAIF*(SIN(BETA)-SIN(BETA+MU) ) 
l+0.5*( (SIN(BETA+MU) )**2-(SIN(BETA) )**2)*PBIF 
FF (2f 2 )=PABETA*(SIN( BETA )-SIN( BETA+MU) ) 

1+A*( COS (BETA) -COS (BETA+MU) ) 

H0.5*( (SIN (BETA+MU) )**2-(SIN(BETA) )**2)*PBBETA 
1 +B* ( BIN ( BETA+MU ) *COS ( BETA+MU) -S1N( BETA ) *COS ( BETA ) ) 

1 +0 . 25*PCBETA* ( 2*MU+SIN ( 2* ( BETA+MU ) ) -S IN ( 2*BETA ) ) 
1+0.5*C*(C0S(2*(BETA+MU) )-C0S(2*BETA) ) 

FF (2f 3 )=PAMU*(SIN( BETA) -SIN( BETA+MU > ) 
1-A*C0S(BETA+MU)+0.5*PBMU*( (SIN(BETA+MU) )**2 
1-(SIN(BETA) )**2)+B*(SIN(BETA+MU)*C0S(BETA+MU) ) 
1+0.25*PCMU*(2*MU+SIN(2*(BETA+MU) )-SIN(2*BETA) > 

1+0,5*C*(1+COS( 2* (BETA+MU) ) ) 

FF(2f 4 )=DELTAO+PAV*(SIN (BETA) -SIN (BETA+MU) ) 
1+0.25*PCV*(2*MU+SIN(2*(BETA+MU) >-SIN( 2*BETA ) ) 

FF ( 2 F 5 ) =PAU* ( SIN ( BETA ) -SIN( BETA+MU) ) 

1+0.5*PBW*( (SIN(BETA+MU) ) **2- ( SIN( BETA ) ) **2 ) 

FF(3f1 )*(3.*0MEGA/3. 1416)*(1 .732*MF*SIN(BETA+MU) ) 
FF(3f2)=(3.*0MEGA/3. 1416)*( ( 1 ,732*IF*HF*C0S( BETA+MU) ) 
l+( (9. *M0/3. 1416)* (W*COS (BETA+MU) -V*SIN(BETA+MU) ) ) ) 
f F(3f3)  = (3.*0MEGA/3. 1416)*( (1 ,732»IF*MF*C0S(BETA+MU) ) 

1 + ( ( 9 . *M0/3 . 1 41 6 ) » ( U*COS ( BETA+MU ) -U»8IN ( BETA+HU ) ) ) ) 

FF ( 3 F 4 ) -27 . *OMEOA«MO»COS ( BETA+MU )/((3.1416)««2) 
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FF  (3»5)=27.*0MEGAlcM0*SIN<BtTA+Mll)/(  (3. 1416)»*2) 

FF (4, 1 )=MF*C0S(BETA)/1 .732 
FF  <4.2)  = (-IF*MF*SIN(BETA)/l .732) 
l+( (3.*M0/3. 1416)* (-W»SIN( BETA )-V*COS( BETA) ) ) 

FF (4.3) =-3.*M0»IL*SIN(MU)/3. 1416 
FF(4.4)=-3.*M0*SIN(BETA)/3. 1416 
FF(4.5)=3.*M0*C0S(BETA)/3. 1416 
FF<5. 1 )=-4.*MF*C0S(BETA+MU/2)*SIN<MU/2)/l .732 
FF(5.2)=(4.*IF*MF*SIN(BETA+MU/2)»SIN<MU/2)/l .732) 
l+( (6.*MO/3.1416)*(W*(COS(BETA)-COS(BETAFMU) ) 
1-V*(SIN<BETA)-SIN<BETA+MU) ) ) ) 

FF(5.3)=(-4.*IF*MF/1 . 732 )*< -0 . 5*SI N ( BF TA+MU/2 ) *S 1 N ( MU/2 ) 
l+0.5*C0S(BETA+MU/2)*C0S(MU/2) ) 

l+(6.*M0/3. 1416)*<-W*C0S(BETA+MU)+V*SIN(BETA+MU) ) 
l-6.*IL*M0*C0S(MU)/3. 1416 
FF(5.4)=6.*M0*(CC)S(BETA)-C0S(BETA+MU)  )/3. 1416 
FF(5.5)=6.*M0*(SIN(BETA)-SIN(BETA+MU) )/3. 1416 
DO  2 11=1.5 


FF(  1 

.11) 

=FF( 

1 . 

I 

I 

)*1 

.OE 

04 

CONT 

INUE 

DO  3 

11  = 

1.5 

FF(2 

.11) 

=FK( 

2. 

I 

I 

)*1 

,0E 

04 

CONT 

INUE 

DO  4 

11  = 

1.5 

FF(4 

. II  ) 

=FF( 

4. 

I 

I 

)»1 

.OE 

04 

CONT 

INUE 

DO  5 

11  = 

1.5 

FF<5 

.11  ) 

=FF( 

5. 

I 

I 

)«1 

.OE 

04 

CONTINUE 

RETURN 

END 


141. 


SUBROUTINE  JAC0B4 ( MF » BETA » IF » W f V» IL »MO»MUrDELTAO»OMEGA»FF» 
1A>B>C) 

DIMENSION  FF(4»4) 

REAl.*4  MOrMOOf  MUf  MFf  ILf  IFfKl  f OMEGA 
PAIF=1 . 155*MF»SIN(BETA) 

PBIF=1 . 155*MF 

PABETA=1 . 155*IF*MF*C0S(BETA)+1 .91*M0*(W*C0S(BETA)-U*SIN(BETA) ) 
PBBETA=-1 .91»M0»IL*SIN(BETA+MU) 

PCBETA=-1 ,91»M0*IL*C0S(BETA+MU) 

PAMU=-1 .91*M0*IL»C0S(MU) 

PWMU=0.0 

PUMU=0.0 

PBMU=~1 .91*M0»IL*SIN(BETA+MU) 

PCMU=-1 .91*M0»IL*C0S(BETA+MU) 

PAW=1 .91*M0*SIN(BETA) 

PBW=1 .91*M0 

PAU=1 .91*M0»C0S(BETA) 

PCU=1 .91*M0 

FF( 1 f 1 )=(COS(BETA+MU)-COS(BETA) )*PABETA 
1-A*(SIN<BETA+MU)-SIN(BETA) )-0,25»(2.*MU 
1-SIN(2.*(BETA+MU) )+SIN(2.*BETA) )»PBBETA 
1+0.54B»(-C0S(2.*(BETA+MU) ) +COS ( 2 . *BETA > ) 

H0.5*(  (SIN(BETA+MU)  ) **2- ( SIN<  BETA ) )*»2) 

1 ♦PCBETA+C*( SIN (BETA+MU)*COS(BETA+MU) -SIN ( BETA )*COS( BETA ) ) 

FF< 1 »2)=PAMU*<C0S(BETA+MU)-C0S(BETA) )-A*SIN(BETA+MU) 
1+0.25*<2.*MU-SIN(2.*(BETA+MU) )+SIN(2.*BETA) )» 
lPBMU+(B/2. )*< 1-C0S(2.*(BETA+MU) ) )+0.5*( ( SIN< BETA+MU ) )*»2 
1-(SIN(BETA) >**2)*PCMU+C*(SIN(BETA+MU)*C0S(BETA+MU) ) 

PP< 1 »3>=PAU*(C0S<BETA+MU)-C0S(BETA) )+0.5*( (SIN(BETA+MU) )»»2 
1- (SIN (BETA) )**2)»PCU 

FF( 1 »4)=DELTA0+PAU*(C0S(BETA+MU)-C0S(BETA) ) 

1 +0 . 25* ( 2 . *MU-S I N ( 2 . * < BETA+MU ) ) +SI N ( 2 . *BETA ) ) *PBW 
FF(2f 1 )=PABETA*(SIN(BETA)-SIN(BETA+MU) ) 
1+A*(C0S(BETA)-C0S(BETA+MU) ) 

1 f0.5*( (SIN(BETAFMU) ) **2- ( SIN ( BETA ) )*»2)*PBBETA 
1 +B* ( S I N ( BETA+MU ) »COS ( BETA+MU  > -SI N ( BET A > *COS ( BET A ) ) 

1 +0 . 25*PCBETA» ( 2 . *MU+SIN( 2 . ♦ ( BETA+MU ) ) -SIN( 2 . *BETA ) ) 

1+0.5*C*( COS (2.* (BETA+MU) ) -COS ( 2 . *BETA ) ) 

FF ( 2 » 2 ) =PAMU* ( SIN<  BETA ) -SIN ( BETA+MU ) ) 
1-A*C0S(BETA+MU)+0.5»PBMU*( ( SIN( BETA+MU ) )t*2 
l-( SIN (BETA) ) »*2 ) +B* ( SIN ( BETA+MU >*COS( BETA+MU) ) 
l+0.25»PrMU*(2.*MU+SIN(2.*(BETA+MU) )-SIN(2.*BETA) ) 

1+0.5»C*( 1+C0S( 2.* (BETA+MU) ) ) 

FF(2f 3)=DELTA0+PAU* (SIN (BETA)-SIN (BETA+MU) ) 

1 +0 . 25*PCU» ( 2 . *MU+S IN ( 2 . * ( BETA+MU ) ) -SIN( 2 . *BETA ) ) 

FF(2f 4) =PAW*(SIN( BETA) -SIN (BETA+MU) ) 

1+0.5*PBW*( (SIN (BETA+MU) ) **2- ( SIN( BETA ) )**2) 

FF(3r 1 )=(3.*0MEGA/3. 1416)*( ( 1 . 732* IF*MF*COS ( BETA+MU ) ) 
l+( (9. *M0/3, 1416)* (g*COS( BETA+MU) -V*S1N(BETA+MU) ) ) ) 
FF(3»2)=(3.*0MEGA/3. 1416)*( ( 1 . 732* IF*MF*COS ( BETA+MU ) ) 
l+( (9.4M0/3. 1416)*a.*C0S(BETA+MU)-V*SIN(BETA+MU) ) ) ) 

FF ( 3 f 3 ) =27 . «0MEGA»M0*C0S ( BETA+MU ) / ( ( 3 . 1416 ) **2 ) 

FF ( 3 » 4 ) =27 . *0MEGA*M0*S I N( BETA+MU ) / ( ( 3 . 1 41 6 > **2 ) 

FF(4.1 )=(4.*IF*MF*SIN(BETA+MU/2. )*SIN(MU/2. )/l .732) 
l+( (6.4M0/3. 1416)*(W*(C0S(BETA)-C0S(BETA+MU) ) 
1-U*(SIN(BETA)-SIN(BETA+MU) ) ) ) 

FF(4»2)=(-4.*IF*MF/1 .732)*(-0.5*SIN(BETA+MU/2. )*SIN(MU/2. ) 
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1 f0.5*C0S(  BFTA  + MU/2.  )*C0S(Ml)/2.  > ) 

l + (A.*M0/3. 1 41A  )*(-U*COS  (E<F1  A fHU)+ V<SIN<E<LFA4  MU)  ) 
l-A.*M0*C0S(MU)*IL/3. 1416 
FF(4»3)=A.»M0*(CC)S(PETA)-C0S(BfcTA+MU)  )/3. 1416 
FF(4»4)=6.«M0<<SIN(  f<E  TA)-SIN(HETA  + MU)  )/3. 1416 
DO  2 11=1 r4 

FF<  1 » 1 1 )=FF( 1 f 1 1 )*1 .OE  04 

CONTINUE 

DO  3 1 1 = 1 » 4 

FF<2»  II  )=FF(2»  II  XI  .OE  04 

CONTINUE 

DO  4 II  = 1 r 4 

FF  < 4,  I I )=FF(4»  II  XI  .OE  04 

CONTINUE 

RETURN 

END 
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SUBROUTINE  LI TLI ( DELTAO » HF r IF » BETA » MO » W . 

IMOOfVr IL»DELTAF»MU»AID»AIQ»AIF»AIK»THETA»KOfKFfKD) 
DIMENSION  AIF(60) r AI D( 60 > > AIQ( 60 > r AIK ( 60) f THETA( 60 ) 

RE AL MF  f I F > MO  > MOO  f I L » MU  f MD  > KQ  r KD  r KF 
WRITE(7f 10)KF»KQ»KD 

10  FORMAK'  'f1X»'KF  ='fE15.3f'KQ  -'fE15.3r'KD  -'rE15.3) 
URITE(7» 1 1 ) IL 

URITE(7.12)BETA 
WRITE(7»13)MU 
URITE<7f14)IF 
URITE(7f 15)0 
URI rE( 7f 16)W 

11  FORMAK'  'rlXf'IL  ='fF10.2) 

12  FORMAK'  'fIXf'BETA  ='fF10.3) 

13  FORMAK'  'fIXf'MU  ='fF10.3) 

14  FORMAK'  'flXf'IF  ='fF10.3) 

15  FORMAK'  'rlXf'O  ='fF10.3) 

16  FORMAK'  'flXf'W  ='»F10.3) 

WRITE(7f9) 

9 FORMAT('  ' »9Xf 'THETA' »9X» ' IK' »12Xf ' 10' r 12Xf ' ID' »12X» 'IF' ) 

ANG=BETAfMU-l .0472 
DO  100  L=1f60 
X=«ANG-BETA 
Y=ANG-(BETA+MU) 

Z=ANG-(BETA+1 .0472) 

IF(X>  2f2f3 
IF(Y)  5f2.4 
IF(Z)  2f2f3 
AIK(L)=0.0 
GO  TO  6 

AIK(L)=( 1 ./DELTAO)*< 1 . 155*MF*IF* ( SIN ( BETA ) -SIN( ANG ) ) 

1+1 .91*M0*U*(SIN<BETA)-SIN(ANG) ) +1 , 91 *MOO*V* < COS ( BETA ) 
l-COS(ANG) )-0.955*IL*DELTAF*(SIN(MU)+SIN(ANG-BETA-MU) ) ) 

6 AIQ(l.)=l . 732*Ka*(0.955*(V-IL*SIN<BETA+MU)  ) 

l + ( II.*SIN(ANG+0.524)-AIK(L)*C0S(ANG)  ) ) 

AIEXD^l  .732»KD«(0.955*(U+IL»C0S(BETA+MU)  ) 
l-(  IL«C0S(ANG+0.524)+AIK(L)*SIN(ANG) ) ) 

AIF(L)=AID(L)*KF/KD 

THE TA(L)=ANG» 180. 0/3. 1416 

AIK(L)=-AIK(L) 

ANG=ANG+0. 01745 

URITE(7»200)THETA(L) »AIK(L) »AIO<L) »AID(L)fAIF(L) 

200  FORMAT<'  ' » 1 X r F 1 4 . 2 » 4E1 4 . 3 ) 

100  CONTINUE 
RETURN 
END 
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SUBROUTINE  NEWTON ( MU. BETA» IF » Wr UfFREOr DELI  AO » IL fVL » M00» K1 . 

IMF. MO. OMEGA. DELTAF.ZETA) 

DIMENSION  FD(5) .F(5. 1 ) .FF(5.5) 

REAL*4  MO.MOO.MU. IF. IL.Kl .MF 
C SLIGHT  ERROR  FOR  IL.NE.O  BUT.LT.35 

KAT-O 
K=0 

IF< IL.GE.35.0)  GO  TO  50 
MU=0.0 

BETA=3. 1416/2. 

IF=( 1 ./(MF*1 .732) )*< (4. 17E-04)»UL-0.75*1L»DELTA0) 

U=0.0 

U=IL 

GO  TO  51 

50  X«SQRT< (18*FREQ*DELTA0»IL)/(4*UL+18*FREQ*DELTA0*IL) ) 

CALL  ARCSIN(X) 

MU=2»X 

2ETA=SIN(MU/2)/(MU/2) 

A»-< ( <3. 1416*DELTA0)/(6*( l-COS(MU) ) *M00 ) ) + < 1 -K1 ) *SIN( MU ) 
l+ZETA*Kl»SIN(MU/2) )/( < 1 -K 1 ) *COS < MU ) +ZETA*K 1 »COS ( MU/2 ) ) 

B.=ATAN(A> 

1F(B.GE.0. ) GO  TO  2 
BETA=3.1416+B 
GO  TO  3 
BETA  = B 

IF=(1/(1 .732*MF»SIN(BETA+MU) ) )*< (VL/(6*FREQ) ) - < 0 . 75*IL*DELTA0 ) - ( 
14.5*IL*ZETA*Kl*DELTAF*SIN(MU/2)/3.1416) ) 

U = IL*Kl*(-C0S<BETATMU>+ZETA*C0S(BETA+MU/2> ) 
V=IL*Kl*(SIN<BETA+MU>-ZETA*SIN(BETA+MU/2) ) 

WRITE(7.510) 

10  FORMATC '0' .5X. '♦*»**FRANKLIN  SOLUTION********') 
URITE(7.300)K.BETA.MU. IF.W.V 
1 CONTINUE 

no  70  K= 1.200 
KK=0 

A=1 . 1 55*IF*MF*SIN( BETA ) T1 .91*M0*<U*8IN( BETA) 

1+0*C0S(BETA) )-l .91*M0*IL*8IN(MU) 

B=1 . 155*IF*MFT1 .91*M0*W+1 .91*IL*M0*C0S( BETA+MU) 

C=1 .91*M0*(V-IL*SIN(BETA+MU) ) 

CALL  RHS( IF.MF.BETA.MO.W.V. IL.MU.DELTAO.F.OMEGA.A.B.C.VL) 

C WRITE(7.486)F(1.1).F(2.1).F(3.1).F<4.1 ).F(5.1) 

486  FORMAT (5E20. 7) 

DO  71  L*1.5 
FD(L)=-F(L.l ) 

Y=ABS(FD(L) ) 

IF(Y.GT.0.001 ) KK»1 

71  CONTINUE 
IF(KK.GE.l)  GO  TO  72 
GO  TO  75 

72  CONTINUE 

CALL  JACOB (MF. BETA. IF. U.V. IL.MO.MU.DELTAO.OMEOA.FF.A.B.C > 
IF(KAT.NE.O)  GO  TO  100 
100  KAT=1 
IEQN-5 
IUEC*1 
EPS»0.01 
KlO-0 
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CALL  0ELG(FDrFF»IEQNrIVEC»EP8rK10) 

BETA»BETA+FD(2) 

HU-MU+FD(3) 

IF*IF+FD(1) 

V=V+FD<4) 

W=W+FD(5) 

C URITE(7f300)KfBETAfMUf IFfUfV 

300  FORMAK'  ' » IX » 14 » 3X » ' BETA- ' f E14 .7»  3X » 'HU* ' »E14 . 7r 

13Xr 'IF-' »E14.7»3X» ' W- ' »E14 , 7 . 3X» ' V- ' »E14 ♦ 7 > 

70  CONTINUE 

gRITE(7.78) 

78  FORMAT ('  ' f IX f ' NEWTON-RHAPSON  DOES  NOT  CONVERGE ' ) 

75  CONTINUE 

IF(K,E0.1)  MRITE(7f500) 

500  FORMAT<'  'f1X» 'NEWTON  DID  NOT  ITTERATEr  K-1') 

C WRITE(7»487) 

487  FORMAT < '0' f2Xf '«*»»$«*««FINAL  SOLUTION  *««»»*»»»') 

C URITE<7f300)KFBETAfMUf IFfUrV 

C URITE<7f488) 

488  FORMAT ( '0' »2Xr 'THIS  IS  THE  RHS  VECTOR  FOR  TEST') 

C WRITE(7»489>F(l»l)fF(2»l)fF<3Fl)»F<4»l)rF<Sfl) 

489  FORMAT (3E20. 7) 

RETURN 

END 


146 


W M 


SUBROUTINE  NEW30N< MU » BETA r IF  »W» V»FREQ» PELTAO» IL »VL r MOOfKl . 

1 MF  f MO » OMEGA . DELTAF » ZET A f K ) 

DIMENSION  FD(5) .F(5» 1 ) ,FF(5»5) 

REAL*4  MOfMOOfMU» IF» ILfKl »MF 
C SLIGHT  ERROR  FOR  IL.NE.O  BUT.LT.35 

KAT=0 

IF< IL.GE.35.0)  GO  TO  50 
MU=0.0 

BETA=3. 1416/2. 

IF=( 1 ./(MF*1 .732) )*( (4. 1 7E -04 ) *UL-0 . 75* I L4DELT AO ) 

U=0.0 

U=IL 

GO  TO  51 

50  X=( (4*3. 1416*VL-9*0MEGA*IL*DELTA0)/<4*3. 1416*VL 
1+9*0MEGA*IL*DELTA0) ) 

CALL  ARCCOS(X) 

MU=X 

ZETA=SIN(MU/2)/(MU/2) 

A=-( ( (3. 1416*DELTA0)/(6*( l-COS(MU) )*M00) ) + < 1 -K 1 ) *SIN ( MU ) 
l+ZETA*Kl*SIN(MU/2) )/( ( 1-Kl ) *COS ( MU > +ZETA*K 1 *COS ( MU/2 ) > 

B=ATAN(A> 

IF<B.GE.O. ) GO  TO  2 
BETA*3. 1416+P 
GO  TO  3 
BETA  = B 

IF=( l/( 1 .732*MF*SIN(BETA+MU) ) )*< ( UL/ <6*FREQ ) ) - ( 0 . 75*IL*DELTA0 ) - ( 
14 .5*IL*ZETA*Kl*DELTAF*SIN(MU/2)/3. 1416) ) 

U = IL*Kl*(-C0S(BETA+MU)+ZETA*C0S(BETA+MU/2) ) 
U=IL*Kl*(SIN(BETA+MU)-ZETA*SIN(BETA+MU/2) ) 

51  CONTINUE 

DO  70  K=l»200 
KK=0 

A=1 . 155*IF*MF*SIN(BETA)+1 . 91 *M0* ( W*SIN( BETA) 

1+V*C0S(BETA) )-l .91*M0*IL*SIN(MU) 

B=1 . 155*IF*MF+1 .91*M0*W+1 .91*IL*M0*C0S(BETA+MU) 

C=1 .91*M0*(V-IL*SIN(BETA+MU) ) 

CALL  RHS( IF»MFfBETA»MO»W»U» IL»MU»DELTAO»FfOMEGA»A»B»C»VL) 

C gRITE(7»486)F( 1 »l)»F(2.1).F(3»l)fF(4fl)FF(5»l) 

486  FORMAT (5E20. 7) 

DO  71  L«l»5 
FD(L)=-F(L» 1 ) 

Y=ABS(FD(L) ) 

IF(Y.GT.0.01 ) KK»1 

71  CONTINUE 
IF(KK.GE.l)  GO  TO  72 
GO  TO  75 

72  CONTINUE 

CALL  JACOB < MF » BETA  r IF  f W f r I L » MO . MU f DELTAO » OMEGA  f FF  f A » B » C ) 

IEQN=5 

IVEC=1 

EPS*0.01 

KlO-0 

CALL  0EL0(FDfFF»IEQN.IVEC»EPS»K10) 

BETA-BETA+FD(2) 

MU-MUfF0(3) 

IF«IF+FD<1) 

V-V+FD(4) 
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U=W+FD(5) 

C URITE<7f300)KrBETA>HUf IFfUfV 

300  FORMAK'  ' f IX » 14 » 3X  f ' BETA- ' f E14 . 7 . 3X » ' MU« ' »E14 . 7r 

13Xf 'IF=' fE14,7f3Xf ' W= ' f E14 . 7r 3X f 'V-' fE14.7> 

70  CONTINUE 

URITE(7f78) 

78  FORMAT<'  ' f IX f ' NEWTON-RHAPSON  DOES  NOT  CONVERGE') 

75  CONTINUE 

IF(K.Ea.l)  WRITE(7f500) 

500  FORMAT <'  'f IX f 'NEWTON  DID  NOT  ITTERATEf  K-1 ' ) 

C WRITE(7f487) 

487  FORMAT< '0' f2Xf '«»««««»*»FINAL  SOLUTION  »»««»«»»»') 

C URITE(7f300)KfBETAfMUfIFfUfV 

C WRITE(7f488) 

488  FORMAT ( '0' f2Xf 'THIS  IS  THE  RHS  VECTOR  FOR  TEST') 

C WRITE <7f 489 )F(1f1)fF(2f1)fF(3f1)fF(4f1)fF(5f1) 

489  FORMAT (5E20. 7) 

RETURN 

END 
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SUBROUTINE  PHACON( IF » MT  » BETAf MO » W » V , IL , MU » DEL T AO r UMEGA » UL » MOO ) 
REAL*4  IFrMFfMOt IL.MU»MOO»OMEGA 
DIMENSION  FDl (4) »F1 (4» 1 ) .FFl (4»4) 

KAT=0 

BETA=BETA+( 5.0*3. 1416/180.0) 

DO  100  K«1f70 
KK=0 

A=1 . 155* IF»MF*SIN( BETA) +1 .91 *MO*(U*SIN( BETA) 

1+U*C0S(BETA) )-l .91*M0*IL*SIN(MU) 

B=1 . 155*IF*MF+1 .91*M0*M+1 .91*IL*M0*C0S(BETA+MU) 

C=1 .91*M0*(0-IL*SIN<BETA+MU) ) 

CALL  RHS4B4< IFfMFfBETAfMOfWrVf ILfMUfDELTAOfFI »OMEGA»A»BrCfUL) 
no  101  L=l»4 
FDl (L)=-F1 (Lf  1 ) 

Y=ABS(FD1 <L) ) 

IF( Y.OE.0.01 ) KK=1 

101  CONTINUE 
IF(KK.GE.l)  GO  TO  102 
GO  TO  105 

102  CONTINUE 

CALL  JAC0B4(MF»BETA» IFfWrO. IL » MO » MU . DEL r AO . OMEGA » FF 1 » Af B . C ) 
IF(KAT.NE.O)  GO  TO  300 
300  KAT=1 
IEQN=4 
IVEC=1 
EPS=0.01 
K10=0 

CALL  GELG(FD1 »FF1 f IEQNf IUECfEPSfKIO) 

BETA-BETA+FDl ( 1 ) 

MU*MU+FD1<2) 

0=V+FD1<3) 

U=-UFFD1  (4) 

100  CONTINUE 

WRITE(7f2000) 

2000  FORMAT('  ' f 1 X f ' NEWT-RAP  DOESN  T CONV  FOR  PHASE  CONTROL') 

105  IF<K.EQ.l)  WRITE(7f305) 

305  FORMAT ('  'fIXf'PHACON  DID  NOT  ITTERATEf  K-1 ' ) 

RETURN 

END 
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SUBROUTINE  RHS( IF rHF » BETArHOf U»Ur IL »MUr DELTAOrFf OHEOAf Ar Br C»UL ) 
DIMENSION  F(5»l) 

REAL«4  MO»HUfMF»IL>IFfOHEOA 

F < 1 » 1 ) =DELTAO»W+A*  < COS  ( BETA^^MU ) -COS  ( BETA ) ) + < B/4 . ) * 

1 (2.«MU-SIN(2.«<BETAFMU)  )-t-SIN<2.«BETA>  ) 
l+(C/2. )»( (SIN(BETA+MU) )**2-(SIN(BETA) )**2> 

F(2f  1 )=DELTAO*U-»^A» (SIN (BETA) -SIN(BETA-I-MU) ) 
lT(B/2. )*( (SIN(BETAFMU) )««2-( SIN (BETA) )»«2) 
lT(C/4. )*(2.*MU+SIN(2.*(BETA+HU))-SIN(2.*BETA) ) 

F ( 3 » 1 ) =- VL+0 . 955*0MEGA* ( 0 . 75* ILBDELTAO+ 1 . 732* IF*HF* 
lSIN(BETAfMU)F2.e65*H0*(U*SIN(BETA-l-MU) 
l+U*COS(BETA+MU) ) ) 

F(4» 1 )=(0,5774*IF)*MF*C0S(BETA)+(0.955*M0) 

1 * ( U«COS  ( BETA ) -U*SI N ( BETA ) ) -t-0 . 9S5*HO*COS ( MU ) * IL 
F ( 5 » 1 ) =-DELTAO*IL- ( 2 . 3094*IF ) *MF*COS ( BETA+MU/2 ) 
l*SIN(MU/2)+( 1 .91*M0)*(W*(SIN(BETA) 

1-SIN(BETA+MU) ) +V* ( COS ( BETA ) -COS ( BETA+MU ) ) ) 

1-1 .91«M0*SIN(MU)*IL 
F(l»l  )=F(1»1)*1.0E  04 
F (2» 1 )=F(2f 1 )*1 .OE  04 
F(4»1)=F(4»1)*1.0E  04 
F(5» 1 )=F(5» 1 )*1 .OE  04 
RETURN 
END 
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SUBROUTINE  RHS4B4 ( IF f MF » BE  I A pMO » W » U r II  f MU» HEL FAO »F  rOMEGAr 
lAfBfCfUL) 

DIMENSION  F(4rl  ) 

REAL*4  MO»MUrMF»IL»IFf omega 

F( 1 »1 )=DELTAO*U+A*(COS(BETAFMU>-COS(BETA) ) + (B/4.  )* 

1 (2.*MU-SIN(2.*(BETA+MU) ) f S I N ( 2 . *BE TA ) ) 
l+(C/2. )»( (SIN(PETA+MU) )**2-( SIN (BETA) )**2) 

F(2. 1 )=DELTAO*U+A»(SIN(BETA)-SIN(BETA+MU) > 
l+(B/2. )»( <SIN(BETA+MU) ) *»2- < SIN ( BETA ) )»»2) 
l + (C/4.  )*(2.*MU  + SIN<2.4(BETA+MU) )-SIN(2.*BETA) > 

F(3»l )=-UL+0.955*0MEGA»(0.75»IL*DELTA0+l .732*IF«Mf  * 
1SIN(BETA+MU)+2.865»M0*(W*SIN(BETA+MU) 
l+U*COS(BETA+MU) ) ) 

F<4fl )=-DELTA0*IL-(2.309*IF)*MF*C0S(PFTA+MU/2. ) 
l«SIN(MU/2. )+( 1 .91*M0)»(W*(SIN(BErA) 

1-SIN(BETA+MU) ) ( COS ( BETA ) -COS ( BETA+MU ) ) ) 

1-1 .91*M0*SIN(MU)*II 
F( 1 » 1 )=F( 1 » 1 )*1 .OE  04 
F(2fl)»F<2»l)*1.0E  04 
F(4»1)=F(4»1)*1.0E  04 
RETURN 
END 


SUBROUTINE  RHS< BETArLLL r MUtRNSIF » IL »DELTAOf MF f MOf Ur V»NOOf IF 
l»KFrDELTAF) 

DIHENSION  RHSIF(50) 

REAL»4  IKflFf IL>MUfMFfHOOtMO»KF 
THE7A=BETA+HU-(3. 1416/3. ) 

DTHETA=*3. 1416/300. 

SIF-0.0 

DO  1 K=lf50 

THETA=THETA+DTHETA 

IK=( 1 ./DELTAO)*<  <2./ 1.732 )*IF*MF*< SIN( BETA) -SIN (THETA > ) 
lT<6./3. 1416)* (HO»U«(SIN< BETA) -8IN< THETA) )iH00*0»( COS (BETA) 

1 -COS ( THETA ) ) )-(3./3.1416)*IL*DELTAF*(SlN(MU)+SIN(THETA- 
IBETA-HU) ) ) 

IF(THETA.LT.BETA)  IK=0.0 

SIF=(1 .732*KF*( (3./3.1416)*(W+IL*C0S(BETA+MU) )-(IL*COS( 
lTHETA+3. 1416/6. )+IK*SIN( THETA ) ) ) )**2+SlF 
THETA=THETA+DTHETA 
CONTINUE 

RMSIF(LLL)»SQRT(SIF)/SQRT(50. ) 

RETURN 

END 
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SUBROUTINE  TERMA ( N» Y »PHI » X » AN! r BNl ) 

AN1-Y*(SIN< (N+1 >»X)*COS(PMI )+COS( ( N+ 1 ) *X ) »S IN ( PHI ) ) 
1/(2*(N+1 ) ) 

ANl>ANl-fY*(SIN<  <N-1  ) *X  ) *COS  ( PHI ) -COS  ( < N- 1 ) *X  ) *SIN ( PHI  ) ) 
« 1/(2*(N-1>) 

BNl=Y*(-COS( (N+1 )»X)*COS(PHI )+SIN( ( N+ 1 ) *X ) »S IN ( PHI ) ) 
1/(2*(N+1 ) ) 

BN1=BN1-Y*(C0S( (N-1 )»X)*COS(PHI )+SIN( (N-1 )»X)*SIN(PHI ) ) 
1/(2*(N-1 ) ) 

RETURN 

END 


t 
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